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14.  ABSTRACT 

The  purpose  of  this  research  project  is  to  develop  new  approaches  that  can  impact  early  diagnosis  of  breast  cancer,  post-biopsy  analysis,  lymph  node  examinations, 
and  drug  delivery  studies.  The  scope  of  this  training  involves  optimizing  Raman  instrumentation  and  methods  for  efficient  illumination  and  collection  of  Raman 
scattered  light  originating  from  deep  within  breast  tissue.  Experiments  encompass  designing  tissue  phantoms  using  Intralipid®,  dyes/pigments,  inclusions,  and 
agarose  gel  to  quantitatively  characterize  the  instruments’  capabilities  and  performance.  Additionally,  we  conduct  spectroscopy  on  tissue  biopsies  to  correlate 
spectral  bands  with  healthy  and  diseased  breast  tissue.  Observed  spectral  changes  are  compared  to  infrared  images  and  H&E/HIS-stained  images  to  relate 
observed  biochemical  changes  to  histology.  The  Raman  data  will  be  used  to  generate  and  train  classification  algorithms  for  automated  histopathology  as  a  starting 
point  for  in  vivo  work  at  the  culmination  of  this  training  program.  This  report  outlines  progress  for  year  1  of  the  3  year  training  program.  In  the  approved  statement  of 
work  5  tasks  were  stated.  All  5  tasks  were  successfully  completed.  We  now  have  a  working  relationship  with  local  clinicians.  We  have  built-up  a  comprehensive 
database  consisting  of  micrographs  and  IR  spectral  images  for  identify  breast  tissue  histology  and  tissue  chemistry.  We  have  acquired  Raman  measurements  on 
biopsied  tissue.  We  have  identified  Raman  spectral  bands  that  can  be  used  for  distinguishing  between  different  tissue  types  and  have  applied  those  Raman 
spectral  bands  to  achieve  cell-level  contrast  in  Raman  spectral  images.  Additionally,  we  have  used  clinical  observations  and  interaction  with  clinicians  to  develop  a 
conceptual  design  for  a  Raman  Tomography  instrument  aimed  at  breast  cancer  screening.  Finally  for  Raman  tomographic  reconstruction,  we  have  evaluated 
existing  diffuse  optical  tomography  algorithms  specifically  for  Raman  measurements  and  have  adapted  a  Monte  Carlo  framework  for  our  Raman  tomographic 
reconstruction.  Through  modeling  and  experimentation  we  characterized  different  instrument  configurations  and  have  implemented  the  transmission  fan  style 
configuration  into  our  prototype  instrument  design  which  uses  fiber-optics  in-contact  with  tissue  to  give  the  maximum  collection  efficiency. 

15.  SUBJECT  TERMS 

Raman  Spectroscopy,  Breast  Cancer,  SORS,  Transmission  spectroscopy,  Raman  Tomography, 

Raman  Imaging,  Infrared  Imaging 
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INTRODUCTION:  The  subject  of  this  research  is  to  develop  non-invasive  optical  instrumentation  and 
methods  for  obtaining  chemical  contrast  from  deep  within  breast  tissue.  The  overarching  goal  of  this 
training  program  is  to  become  steeped  in  breast  cancer  clinical  practice  and  research  to 
independently  develop  a  new  technology  from  a  concept.  My  learning  objectives  include  the 
following:  first,  fill  in  the  gaps  in  my  basic  knowledge  of  breast  cancer  and  pathology.  Second,  build  a 
knowledge  base  as  to  the  chemical  changes  that  are  occurring  within  breast  tissue  and  how  those 
changes  are  reflected  in  the  vibrational  spectra.  Third,  to  understand  the  optics  of  breast  tissue  to 
gain  most  efficient  illumination  and  collection  schemes  for  instrumental  design.  The  purpose  of  this 
research  project  is  to  develop  new  approaches  that  can  impact  early  diagnosis  of  breast  cancer,  post¬ 
biopsy  analysis,  lymph  node  examinations,  and  drug  delivery  studies.  Additionally,  this  technology 
could  be  used  to  study  animal  models  for  tissue  engineering,  3D  tissue  scaffolds  and  has  a  variety  of 
other  basic  research  as  well  as  clinical  applications.  The  scope  of  this  training  involves  optimizing 
Raman  instrumentation  and  methods  for  efficient  illumination  and  collection  of  Raman  scattered 
light  originating  from  deep  within  breast  tissue.  Experiments  encompass  designing  tissue  phantoms 
using  Intralipid®,  dyes/pigments,  inclusions,  and  agarose  gel  to  quantitatively  characterize  the 
instruments'  capabilities  and  performance.  Additionally,  we  conduct  spectroscopy  on  tissue  biopsies 
to  correlate  spectral  bands  with  healthy  and  diseased  breast  tissue.  Observed  spectral  changes  are 
compared  to  infrared  images  and  H&E/HIS-stained  images  to  relate  observed  biochemical  changes  to 
histology.  The  Raman  data  will  be  used  to  generate  and  train  classification  algorithms  for  automated 
histopathology  as  a  starting  point  for  in  vivo  work  at  the  culmination  of  this  training  program. 
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BODY:  The  approved  statement  of  work  is  illustrated  in  figure  1.  For  year  one,  five  tasks  were 
outlined.  All  tasks  have  been  completed  and  we  are  able  to  begin  the  work  outlined  for  year  2  as 
scheduled.  Details  of  the  work  conducted  in  year  one  are  present  on  the  following  pages. 


Task  1 .1  Obtain  tissue  biopsy  slides  and  put  together  a  measurement  plan  (month  1  -2) 

Task  1 .2  Raman  measurements  (month  3-1 0)  |vf 

Task  1 .3  Raman  data  processing  (month  11-12) 


Classification  Algorithm  (month  19-30) 


Year  1 


Year  2 


Year  3 


o 

rsl 


a 

< 


>  Task  2.5  Modifications  to  improve  the 
instruments  performance  (month  25-30) 
►  Task  2.4  Collect  data  using  more  realistic  tissue 
phantoms  (month  19-24) 


b^-Task  2.3  Construction  of  Instrumentation  and  calibration  standards  (month  13-18) 

-►Task  2.2  Instrument  optical  design  in  Zemax©/Matlab®  (month  10-12)  [?f 
Task  2.1  Instrument  conceptual  design  based  on  clinical  observation  (month  3-12) 

Figure  1)  Summary  figure  and  timeline  for  approved  statement  of  work 


Task  1.1  Obtain  tissue  biopsy  slides  and  put  together  a  measurement  plan  (month  1-2):  This 
task  was  successfully  completed.  Tissue  samples  were  obtained  from  Provena  covenant  medical 
center  and  US  Biomax,  Inc  (IRB  approved  sources).  101  breast  tissue  biopsies  cores  (1mm)  were 
obtained.  Diagnosis  were  hyperplasia  (n=20),  Dysplasia  (n=20),  Melignant  Tumor  (n=40),  Normal 
(n=ll).  Serial  sections  of  these  biopsies  were  stained  and  imaged  to  establish  a  database  of 
micrographs  that  could  be  used  to  identify  cell  types  and  regions  of  interest  in  the  tissue.  The  stains 
used  for  identifying  tissue  histology  were  carefully  chosen  upon  consultation  with  two  pathologists. 
The  histological  stains  used  were  as  follows: 


No  Stain  For  Raman  -Serial  Section  261 
P53-Serial  Section  262 
Ki67-Serial  Section  263 
CD31-Serial  Section  264 
P63-Serial  Section  265 
No  Stain  For  IR-Serial  Section  266 
H&E-Serial  Section  267 
Calponin  -Serial  Section  268 


Masson's  Trichrome-Serial  Section  269 
Vimentin-Serial  Section  270 
Smooth  Muscle  Actin-Serial  Section  271 
HMW  Cyto keratin-Serial  Section  272 
No  Stain  For  IR-Serial  Section  273 
Her2/neu-Serial  Section  274 
Estrogen  Receptor-Serial  Section  275 
Progesterone  Receptor-Serial  Section  276 
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For  each  biopsy  we  built  a  database  of  micrographs  that  can  be  zoomed  to  a  20x  magnification.  A  set 
low  magnification  serial  sections  are  depicted  in  figure2  showing  the  different  staining  results. 


HMWCK  Her2/neu  ER  PR 


Figure  2)  Serial  sections  and  staining  obtained  for  each  breast  tissue  biopsy 


In  addition  to  stained  serial  sections  we  also  acquired  mid-infrared  images  from  each  of  these 
biopsies.  Some  of  our  mid  infrared  results  is  reported  in  the  manuscript  'High  Definition 
Spectroscopic  Imaging'  attached  in  the  appendix  of  this  report.  An  example  of  an  infrared  spectral 
image  is  depicted  in  figure  3.  These  images  are  useful  in  working  with  pathologists  in  order  to 
evaluate  histology  and  pathology.  We  used  this  electronic  database  of  infrared  and  white  light 
micrographs  to  identified  regions  of  interest  for  Raman  measurements.  Regions  of  interest  were 
selected  by  correlating  the  tissues  chemistry  with  the  observed  tissue  staining  to  select  different 
tissue  types  i.e.)  distinguishing  epithelial  cells  from  myoepithelial  cells. 


Mid-IR  Tissue  Spectrum 


Figure  3)  a)  Infrared  image  of  using  the  Amide  A  band  for  contrast  b)  representative  mid  IR  spectrum 
illustrating  the  information  available  for  contrast  at  each  pixel  in  this  image. 
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Task  1.2  Raman  measurements  (month  3-10):  We  have  successfully  completed  this  task.  Raman 
spectra  were  acquired  on  a  Horiba  LabRam  microscope  using  a  video  reference  to  mark  the  position 
of  each  measurement.  Spectra  were  obtained  using  785nm  excitation  with  a  lOOx  0.90NA 
microscope  objective.  Acquisitions  range  from  15  to  60  seconds  per  measurement  point  depending 
on  signal  to  noise.  Spectra  were  acquired  at  a  spectral  resolution  between  5  and  10  cm1. 


Raman  Shift  (cm'1) 

Figure  4)  Raman  spectra  were  collected  from  select  regions  within  each  biopsy  core  and  compared  based  on 

tissue  type. 


Raman  Band 
(cm1) 


781 


828 


854 


1032 


1085 


1126 


1156 


1310 


1333 


Major  Assignments 


C-C  twisting  mode  phenylalanine 


C-C  twisting  mode  tyrosine 


Cytosine/Uracil  ring  breathing 


Ring  breathing  tyrosine/O-P-O  stretch 


Tyrosine  ring  breathing 


C-C  stretch  Proline,  Valine,  and  protein 
backbone  (a-helix) 


Symmetric  ring  breathing  mode  phenylalanine 


C-H  in  plane  bending  mode  phenylalanine 


C-N  stretch 


C-N  stretch 


C-C  &  C-N  stretching 


Tryptophan  &  Phenylalanine  n(C-C6H5) 


Amide  III 


CH3CH2  twisting  mode 


CH3CH2  wagging  mode 


CH2  bending  mode 


Weighted  towards 

Epithelium 

Epithelium 

Epithelium 

Epithelium 

Stroma 

Epithelium 

Epithelium 

Epithelium 

Stroma 

Epithelium 

Epithelium 

Epithelium 

Stroma 

Epithelium 

Epithelium 

Stroma 


core  is  depicted  in  figure  4.  All  spectra  in  figure  4  were 
phenylalanine  band  at  1003  cm1. 


We  acquired  Raman  spectra  from 
strategic  locations  on  all  biopsy 
cores  that  were  available  on 
serial  section  261  (see  page  5). 
We  also  acquired  Raman 
measurements  from  surgical 
resections  that  were  provided 
through  Provena  medical  center. 
The  primary  focus  of  our 
measurements  was  in 
distinguishing  different  tissue 
types  in  anticipation  of 
developing  an  algorithm 
classifier.  Figure  4  illustrates  an 
example  of  point  measurements 
collected  on  regions  of  interest 
within  a  biopsy  section.  A 
representative  example  of  the 
data  collected  for  each  biopsy 
baselined  and  normalized  to  the 


Figure  4  illustrates  very  clearly  that  different  tissue  types  can  be  identified  by  observing  the 
differences  in  Raman  bands.  For  example,  the  black  spectrum,  collected  over  stromal  tissue,  has  a 
much  larger  amide  III  Raman  band  (1244  cm1  and  1274  cm1)  than  epithelial  or  myoepithelial  tissue. 
This  band  is  correlated  with  a  C-N  and  N-H  vibrational  modes  and  can  be  used  for  contrast  to 
distinguish  epithelial  tissue  from  stromal  tissue.  Additionally  the  band  at  1045  cm1  also  has  a  larger 
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relative  contribution  in  the  stromal  and  myoepithelium  tissue  over  epithelium.  Some  specific 
spectral  band  weightings  are  listed  in  the  table  below  figure  4.  Validating  the  exact  position  of  a 
measurement  was  difficult  as  the  tissue  structure  changed  slightly  in  sequential  serial  sections  and  as 
a  result  this  contrast  is  some-what  difficult  to  visualize  by  point  measurements. 

To  overcome  this  difficulty,  we  have  also  begun  acquiring  Raman  spectral  images  that  allow  us  to 
more  reliably  identify  where  a  Raman  spectrum  originates  from  relative  to  the  surrounding  tissue. 
Using  the  Raman  spectrum  for  contrast  and  having  the  positional  correspondence  of  many  Raman 
spectra,  we  are  able  to  visualize  the  tissue  structure  and  pull  out  average  differences  in  tissue  types. 
While  we  have  begun  work  on  a  tissue  classifier,  we  anticipate  continued  effort  in  the  Raman  imaging 
approach  to  be  useful  in  both  building  and  validating  our  classification  algorithms. 

An  example  of  a  Raman  spectral  image  is  illustrated  in  figure  5.  The  top  right  image  illustrates  a 
Raman  spectral  image  over-laid  on  an  unstained  biopsy  micrograph.  The  primary  contrast  in  this 
image  is  achieved  using  the  band  at  1045  cm1.  By  correlating  the  position  of  the  image  with  stained 
biopsy  serial  sections  in  our  database,  it  is  clear  that  we  can  distinguish  epithelium,  stroma,  and 
myoepithelial  cells  without  dyes  or  label.  This  is  a  significant  finding  because  myoepithelial  cells 
makeup  the  basal  layer  of  normal  mammary  epithelial  tissue.  Their  identification  has  particular  diagnostic 
value  as  they  are  lost  in  malignancy  but  retained  in  most  benign  lesions.  Quantification  of  these  cells  in  vivo 
could  play  a  key  role  in  identifying  invasive  carcinomas. 


Figure  5)  Raman  spectral  image  compared  to  stained  biopsy  sections 

Task  1.3  Raman  data  processing  (month  11-12):  Spectral  processing  involved  calibrating  the 
wavelength  axis  against  a  neon  atomic  emission  lamp  and  using  a  NIST  traceable  white  light  source 
to  correct  for  the  instruments  wavelength  response.  The  Raman  spectra  were  then  imported  into 
Matlab®  and  examined  for  cosmic  rays  (spikes).  If  spikes  were  found  they  were  removed  by  using  a 
median  filter  on  adjacent  points  to  determine  an  interpolating  reference.  The  spike  were  then 
deleted  and  replaced  by  interpolating  data  points  from  both  sides  of  the  spike  to  fill  in  the  gap. 
Baseline  points  on  the  data  were  identified  and  used  for  all  spectra  analyzed.  Using  these  same 
baseline  points  a  multi-linear  baseline  was  then  fit  to  the  data  and  subtracted  to  provide  baselined 
Raman  spectra.  The  spectra  were  then  normalized  either  by  overall  area  or  to  specific  band  areas. 
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Several  algorithms  and  visualization  approaches  were  developed  for  this  work  and  have  been 
reported  in  the  manuscript  "Designing  Transfer  Functions  for  Interactive  Data  Mining  of 
Hyperspectral  Images"  which  is  attached  in  the  appendix  of  this  report. 

Task  2.1  Instrument  conceptual  design  based  on  clinical  observation  (month  1-2): 

Dr.  Krishna  Tangella  M.D.  who  is  Teaching  Faulty  at  the  University  of  Illinois  at  Urbana-Champaign 
College  of  Medicine  has  spent  the  year  with  me  in  a  one-on-one  course  in  Breast  Cancer  pathology. 
Over  this  year  I  have  built-up  my  fundamental  understanding  of  breast  cancer.  I've  worked  closely 
with  the  pathology  team  at  Provena  Covenant  medical  Center  in  Urbana,  IL.  In  working  with  this 
team  I've  been  fortunate  to  observe  many  aspects  of  clinical  practice  associated  with  breast  cancer. 
I've  observed  screening  with  MRI  and  mammography  as  well  as  ultrasound  guided  biopsies.  I've 
seen  the  processing  that  occurs  when  a  biopsy  is  section  and  prepared  for  histological  staining. 
Additionally,  I've  observed  grossing  procedure  conducted  on  tissue  from  mastectomies.  I  was  also 
fortunate  to  attend  a  week  long  intensive  Breast  Imaging  course  covering  clinical  breast  imaging  that 
was  provided  by  the  International  Institute  of  Continuing  Medical  education  (iiCME).  This  course 
was  extremely  useful  in  understanding  current  radiology  practices  and  was  very  motivating  as 
spectroscopy  was  brought-up  as  a  future  technology  of  interest! 

Based  on  these  observations  and  experiences,  we  developed  conceptual  design  for  a  Raman  based 
screening  instrument  as  depicted  in  figure  6. 


Figure  6)  a)  Photo  taken  at  Provena  Covenant  Medical  Center  of  a  mammogram  instrument  b)  Conceptual 
Design  and  AutoCAD  files  for  a  Breast  Cancer  Screening  Raman  Instrument 


The  instrument  design  is  comprised  of  a  stand  that  couples  two  breast  paddles  to  an  optical  table. 
The  stand  has  vertical  and  rotational  motion  to  adjust  to  patient  height  and  allow  for  multiple 
projections  through  the  breast  tissue.  Both  breast  paddles  are  removable  and  can  be  made 
compatible  with  existing  mammogram  instruments.  The  bottom  breast  paddle  has  a  transparent 
glass  plate  that  is  used  to  compress  the  breast  tissue  and  will  allow  the  Raman  excitation  source  to 
illuminate  the  bottom  of  the  breast.  The  light  transmitted  through  the  breast  is  collected  by  fiber 
optics  on  the  top  paddle  and  is  then  relayed  to  an  imaging  Raman  spectrometer.  This  data  is  then 
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used  reconstruct  the  size,  shape,  and  position  of  features  using  the  innate  chemical  contrast  in  the 
breast  tissue.  For  this  training  program  years  2  and  3,  we  will  build  and  test  this  instrument  using 
tissue  phantoms  (physical  breast  models  that  have  similar  optical  properties  to  breast  tissue). 

Task  2.2  Instrument  optical  design  (month  3-12): 

A  transmission  approach  to  illumination  and  collection  has  proven  to  be  the  most  robust  approach  to 
diffuse  optical  tomography  where  as  other  methods,  a  back  scattered  configuration  for  example,  are 
known  to  have  reconstruction  artifacts  like  higher  surface  weightings  of  reconstructed  targets  W.  To 
explore  and  understand  the  use  of  transmission  Raman  measurements  in  tomographic 
reconstructions  we  conducted  a  series  of  experiments  and  modeling  to  aid  in  instrument  design.  We 
compared  two  forward-modeling  methods,  radiative  transport  calculation  (via  Nirfast,  an  open- 
source  diffuse  optical  tomography  modeling  package)  and  Monte  Carlo  simulation  (written  in-house), 
for  the  modeling  of  light  fluence  in  the  phantom.  Reconstruction  of  the  size  and  position  of  buried 
targets  was  attempted  via  an  iterative  modified-Tikhonov  minimization  algorithm  without  the  use  of 
spatial  priors.  The  results  are  validated  against  computed  tomography  (CT)  images  of  the  same 
samples. 

Phantom  Specimens  and  Fabrication 

Tissue  phantoms  were  fabricated  by  dissolving  1  gram  of  agar  (Sigma-Aldrich)  in  50  ml  of  water  at 
95°C.  Different  volumes  of  20%  Intralipid  (Sigma),  an  oil  emulsion  that  is  commonly  used  to  mimic 
the  scattering  properties  of  fatty  tissue/2' 3* 4]  were  added  to  increase  the  scattering  potential.  0.25ml 
or  1.25ml  per  50ml  of  water  were  added  to  produce  phantoms  of  total  Intralipid  concentration  0.1% 
or  0.5%  by  volume.  A  0.1%  concentration  of  Intralipid  represents  the  lower  limit  of  what  is  opaque 
to  the  naked  eye,  while  a  0.5%  solution  close  to  the  scattering  level  exhibited  by  epithelial  tissue. i5'6] 
The  solution  was  poured  the  wells  of  a  cell  culture  plate.  One  PTFE  sphere  of  1/8"  diameter  or  two 
spheres  of  1/16"  diameter  were  suspended  in  each  well  by  thin  histology  needles,  and  the  plate  was 
allowed  to  cool  in  a  refrigerator  for  approximately  one  hour.  The  resulting  phantoms  were 
cylindrical,  1  cm  in  diameter  and  approximately  2  cm  tall.  Each  sphere  was  positioned  such  that  its 
middle  sat  approximately  1.2  cm  above  the  base  of  the  phantom.  The  radial  position  of  the  spheres, 
as  well  as  the  distance  between  two  spheres  embedded  in  the  same  phantom,  were  experiment- 
dependent. 

Raman  Instrumentation 

The  instrument  comprised  of  a  785  nm  excitation  laser  (harrowed  for  these  experiments)  and  CCD- 
coupled  spectrometer  with  a  pre-stage  notch  filter,  a  fiber-optic  bundle  detector  and  a  motorized 
rotation  stage,  as  shown  in  Figure  7.  In  all  cases,  the  laser  and  detector  were  held  at  a  constant  height 
at  the  level  of  the  PTFE  target  while  the  cylindrical  phantom  was  rotated  about  its  vertical 
(symmetrical)  axis.  This  allowed  for  a  single  source/detector  combination  to  sample  the  phantom 
from  multiple  angles,  a  requirement  of  optical  tomography.  iy] 

The  400mW  CW  laser  at  785nm  (Invictus,  Kaiser)  was  fiber-launched  and  passed  through  a 
collimator  to  produce  a  1.25  mm  diameter  beam.  In  all  cases,  this  beam  was  centered  on  the  waist  of 
the  phantom  at  the  height  of  the  target.  Collection  was  performed  via  two  achromatic  doublets  which 
focused  the  photons  exiting  from  the  phantom  onto  a  50-fiber  bundle  (100  um  diameter  each,  in  an 
approximately  5x10  array)  (Fibertech  Optica,  Kitchener,  Ontario,  Canada).  The  focal  lengths  of  these 
achromats  were  60mm  (focusing)  and  either  75mm  or  200mm  (collection)  depending  on  the 
experiment.  The  output  of  the  fiber  bundle  was  imaged  onto  a  spectral  CCD  (iDUS,  Andor 
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Technology;  Belfast;  Northern  Ireland).  The  motorized  rotation  stage  (NR360S,  ThorLabs,  Newton, 
New  Jersey),  controlled  by  LabView  (National  Instruments,  Austin,  Texas),  was  stepped  clockwise  in 
4.5  degree  increments  in  all  experiments  for  a  total  of  80  discrete  sampling  angles. 


Figure  7)  a.)  and  b.)  Photos  and  c.)  schematic  of  Raman  tomography  instrumentation.  Both  the  illumination 
and  collection  fibers  are  centered  on  the  phantom  at  the  height  of  the  PTFE  sphere  target. 

Experiments 

For  phantoms  containing  a  single  target,  three  separate  instrumentation  configurations  (as  shown  in 
Figure  7  were  considered.  In  the  first,  shown  in  Figure  8a,  the  source  and  detector  were  kept  at  a 
180°  angle  while  the  phantom  was  rotated  through  80  positions  to  demonstrate  the  insufficiency  of 
using  a  single  source-detector  angle.  Acquisition  time  was  three  minutes  per  step  (two  frames  of  90 
seconds  each  to  aid  in  data  processing).  The  second  experiment  involved  varying  the  angle  between 
the  source  and  detector,  shown  in  Figure  8b.  In  addition  to  180  degrees,  the  source  was  moved 
relative  to  the  detector  to  form  angles  of  135,  90,  and  45  degrees.  This  setup  mimics  the  “circular  fan- 
beam"  geometry  which  is  commonly  used  in  diffuse  optical  tomography  experiments!1]  For  each 
source  position,  the  stage  was  rotated  through  the  full  80  steps  and  two  spectra  of  30  seconds  each 
(as  opposed  to  90  seconds)  were  acquired  at  each  step.  In  the  final  configuration,  the 
source/detector  angle  was  returned  to  180  degrees,  but  the  achromatic  collection  lens  nearest  to  the 
phantom  was  replaced  with  one  of  focal  length  200mm  as  shown  in  Figure  8c.  The  entire  detector 
setup  was  moved  an  appropriate  distance  away  from  the  phantom  to  keep  the  fibers  in  focus.  This 
expanded  the  collection  region  from  an  area  of  approximately  0.97  mm  x  1.3  mm  to  2.9  mm  x  3.8  mm 
with  the  intent  of  being  able  to  collect  spatial  information  from  each  individual  fiber.  This  fan-like 
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a) 


Collection  V  y  Illumination 


b) 


Multi-angle 

Illumination 


longer  focal  length 

Figure  8)  Three  instrumentation  configurations 
corresponding  to  the  three  experiments  performed,  a) 
Illumination  and  collection  are  fixed  at  a  180  degree 
angle.  Measurements  are  taken  as  the  phantom  is  rotated 
in  80  steps  of  4.5°  each.  The  fiber  bundle  is  treated  as  a 
point  detector  for  80  total  measurements,  b)  Multiple 
scans  are  taken  with  the  source  forming  four  positions 
with  the  detector  (180°,  135  °,  90  0  and  45°).  The  total 
number  of  measurements  taken  is  320  (4  x  80).  c) 
Illumination  and  collection  are  fixed  at  a  180°  angle  as  in 
setup  a),  but  the  focal  length  of  the  collection  has  been 
increased  from  75  mm  to  200  mm.  As  a  result,  the  total 
area  of  collection  is  larger  and  the  fiber  bundle  can  no 
longer  be  treated  as  a  point  detector  -  the  signals  from 
(ten)  individual  fibers  were  evaluated  for  800  total  data 
points. 


geometry  was  chosen  to  mimic  the  cone- 
beam  transmission  computed  tomography 
(CT)  source-detector  configuration.^] 

A  series  of  experiments  involving 
phantoms  with  two  PTFE  spheres  were 
also  performed.  Using  the  instrument 
setup  with  the  longer  focal  length 
collection  (experiment  seen  in  figure  8c) 
and  an  additional  30  mm  focal  length  lens 
to  focus  the  collimated  illumination  beam 
in  an  approximation  of  a  point  source  (see 
later  discussion),  the  response  from  two 
spheres  was  investigated.  Additionally,  a 
third  series  of  experiments  in  which  the 
Intralipid  content  of  the  two-sphere 
phantoms  was  increased  to  0.5%,  was 
performed  in  the  same  manner  with  two 
spheres  and  a  longer  focal  length.  The 
quality  of  the  reconstructions  was 
evaluated  by  comparison  to  micro-CT 
images  showing  absolute  sphere 
positions. 

Data  Processing 

Zemax  was  used  to  model  some  of  the 
optical  configurartions,  while  matlab  was 
found  to  be  more  useful  in  the 
tomographic  reconstructions.  All  data 
processing  was  done  in  Matlab  R2009b 
(The  MathWorks,  Nantucket  MA).  A 
median  filter  was  applied  to  all  collected 
data  to  correct  for  the  cosmic  ray  'spikes' 
which  appear  on  the  CCD.  The  two  Raman 
bands  of  interest  are  PTFE  (732  cm1)  and 
Intralipid  (-800-820  cm1).  As  with  the 
biopsy  Raman  measurements,  a 
multilinear  “rubber  band"  baselining 
procedure  was  used  to  remove  the 
background,  and  the  area  under  each 
band  was  calculated.  During 

reconstructions,  both  the  absolute 


intensity  of  the  PTFE  signal  and  that  of  the 
Intralipid  signal  are  utilized.  The  intensity  of  the  Intralipid  band  is  directly  proportional  to  the 
number  of  photons  reaching  the  detector  for  that  measurement,  and  can  be  considered  an  internal 
standard  which  accounts  for  any  discrepancies  in  the  distance  between  source  and  detector  or  in  the 
amount  of  light  collected.  In  high-scattering  samples  where  the  Intralipid  band  is  obscured  by  noise, 
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the  band's  value  is  approximated  as  being  proportional  to  the  average  intensity  of  the  entire  CCD 
chip. 

Computerized  Tomography  Reconstruction 

Nirfast  [www.dartmouth.edu/~nir/nirfast],  is  an  open-source  software  package  for  simulation  of 
DOT  experiments  and  for  DOT  reconstructions  from  real  or  simulated  datal9]  Based  on  a  finite- 
element  method  for  calculating  photon  fluence  though  specified  geoHmetries,  the  two  processes  that 
can  be  modeled  (NIR  absorbance  and  fluorescence  emission)  can  be  seen  as  analogous  to  Raman 
scattering  events  (conversion  of  Rayleigh  scattered  to  Stokes  scattered  light).  While  the  physical 
processes  behind  these  two  are  not  the  same,  they  are  similar  enough  to  warrant  the  use  of  Nirfast  as 
a  first  approximation  for  2D  reconstructions  from  experimental  Raman  data.  For  reconstructing 
experimental  data,  a  circular  finite  element  mesh  was  generated  within  the  program.  Source  and 
detector  locations  were  added  to  match  each  experimental  setup,  and  the  experimental  data  were 
ordered  and  normalized  to  match  the  range  of  the  signal  from  a  forward  simulation  of  a  'blank'  mesh 
with  no  absorbing  target. 

Also  evaluated  was  a  Monte  Carlo  simulation  written  in-house  for  modeling  photon  fluence.  Using  a 
2D  pixel  mesh,  the  reconstruction  mathematics  and  assumptions  conserved  between  the  two 
methods.  Although  more  computationally  intensive,  Monte  Carlo  simulations  allow  the  user  to 
overcome  certain  assumptions  made  by  the  radiative  transport  equation  (including  the  use  of  a  non- 
idealized  point  source,  implementation  of  geometrically-accurate  detectors,  and  photon  migration 
that  does  not  follow  the  diffusion  regime)  which  may  not  be  appropriate  for  certain  Raman 
experiments.  Monte  Carlo  studies  which  simulate  the  generation  and  propagation  of  Raman  photons 
through  turbid  media  have  been  performed  recently!10'  12i  In  a  similar  manner,  we  compare  the 
results  of  Monte  Carlo  simulations  to  those  of  Nirfast  to  show  that  the  physical  processes  that  govern 
Raman  tomography  experiments  differ  significantly  from  DOT  experiments  and  that  these  two 
modalities  need  to  be  treated  as  separate  entities. 

Results  from  these  experiments 

A:  Radiative  transport  calculations  fNIRFASTl 

Figure  9a  illustrates  the  simplest  configuration  and  the  recorded  Raman  signal  of  a  PTFE  sphere 
relative  to  Intralipid  as  a  function  of  rotation  (Figure  9b).  Two  maxima  are  observed:  one  when  the 
target  passes  the  source  and  another  when  the  target  passes  the  detector.  These  two  maxima  are 
spaced  by  approximately  180°,  which  is  explained  by  uniform  rotation  of  the  target  through  360°.  In- 
between  the  two  maxima,  no  Raman  signal  is  observed;  the  collected  photon  count  at  these  points  is 
below  the  spectrograph's  limit  of  detection. 

Each  of  the  50  fibers  from  the  detector  is  mapped  onto  a  different  region  of  the  CCD  camera,  and  the 
response  from  each  fiber  can  be  individually  determined.  The  behavior  of  individual  fibers 
(detectors)  set  along  the  waist  of  the  phantom  would  be  similar  to  the  results  seen  in  Figure  9b:  the 
PTFE  target  will  pass  by  each  fiber  and  the  source  during  rotation,  resulting  in  two  maxima  with 
slight  position  shifts  as  a  result  of  fiber  position.  If  the  fibers  are  all  packed  closely  together,  we 
expect  this  difference  to  be  minimal;  if  the  fibers  are  spaced  at  larger  intervals,  this  difference  will  be 
much  more  pronounced.  Figure  9c  plots  each  fiber's  (48  total)  center  of  gravity  for  the  second 
maximum  (when  the  target  is  closest  to  the  detector)  as  a  function  of  phantom  rotation.  The  fibers 
show  a  standard  deviation  of  3.1  degrees,  indicating  that  the  majority  of  the  fibers  'see'  the  target 


13 


□  — 


Figure  9)  a)  Schematic  of  180°  setup,  including  the 
direction  of  rotation  and  the  approximate  starting  position 
of  the  PTFE  target  relative  to  the  source  and  detector,  b) 
PTFE  to  Intralipid  Raman  photon  ratio  as  a  function  of 
phantom  rotation,  average  of  all  fibers.  Distinct  increases 
in  signal  occur  when  the  target  passes  the  source  or  the 
detector,  c)  Center  of  gravity  measurements  for  each 
fiber-optic  in  the  detector  bundle.  Due  to  the  small 
spacing  of  fibers,  little  variation  is  seen  in  the  position  of 
the  second  peak  as  seen  by  each  individual  fiber,  d) 
Reconstruction  from  Nirfast.  Data  predicts  two  targets  at 
opposite  ends  of  the  phantom. 


Figure  10)  a)  Schematic  for  the  multi-angle 
setup,  including  the  direction  of  rotation  and  the 
approximate  starting  position  of  the  PTFE  target 
relative  to  the  sources  and  detector,  b)  PTFE  to 
Intralipid  Raman  photon  ratio  as  a  function  of  both 
phantom  rotation  and  source/detector  angle.  Signal 
maxima  occur  when  the  target  passes  near  the  source 
or  the  detector.  As  expected,  the  four  peaks 
representing  the  target  passing  the  detector  are 
separated  by  about  45  degrees  each.  c) 
Reconstruction  from  Nirfast.  Even  though  there  are 
shadows,  the  data  correctly  predicts  the  location  of 
the  target.  These  shadows  are  due  to  normalization 
errors  and  deviations  from  point  source  illumination. 


pass  by  within  six  degrees  of  rotation.  This  small  variation  is  to  be  expected  because  the  assembly's 
collection  area  spans  one  square  millimeter,  contributing  a  very  small  angular  response  of  less  than 
three  degrees.  The  two  maxima  (target  near  source  and  target  near  detector)  are  similar  in  height, 
width,  and  shape.  Without  angular  information  it  is  not  obvious  how  to  distinguish  one  from  the 
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other.  Performing  a  2D  reconstruction  on  the  data  with  Nirfast,  as  shown  in  Figure  9d,  confirms  this 
observation:  the  program  predicts  two  targets  separated  by  180  degrees.  This  result  reinforces  the 
well-accepted  notion  that  multiple  collection  angles  are  needed  to  accurately  determine  the  position 
of  targets  buried  in  scattering  analyte.  W 

Figure  10  displays  Raman  signal  versus  rotation  angle  for  multiple  source  locations  ("circular  fan- 
beam"  geometry).  The  data  have  been  normalized  to  the  second  maxima  (when  the  phantom  is 
nearest  to  the  detector).  As  in  the  previous  example,  each  full  rotation  shows  two  maxima 
corresponding  to  when  the  target  passes  the  source  and  the  detector.  While  the  180°  configuration 
shows  a  separation  of  roughly  180°,  the  other  angles  show  corresponding  separations  close  to  135°, 
90°  and  45°.  It  is  clear  that  the  maxima  corresponding  to  the  target  passing  the  detector  do  not 
change  in  position.  The  Nirfast  reconstruction  in  Figure  10c  supports  the  observation  that  multiple 
illumination/collection  angles  gives  greater  insight  as  to  the  location  of  the  target. 

a) 

- Simulated  -  point  source 


Figure  11)  a)  Illumination  fluence  at  the  phantom  border  at  four  different  angles  relative  to  the  point  of 
illumination.  Solid  line:  calculated  fluence  from  the  radiative  transport  equation.  Dashed  line:  measured 
fluence  using  a  tissue  phantom  (0.1%  Intralipid,  no  targets),  b)  Illumination  structure  for  an  idealized  point 
source,  c)  Ray  tracing  simulation  showing  the  illumination  structure  for  a  directionalized  laser  source. 


Although  one  target  is  predicted,  'shadowing'  occurs  in  regions  not  occupied  by  the  target.  The 
primary  cause  of  this  stems  from  a  discrepancy  between  how  illumination  is  modeled  and  how  the 
phantom  is  actually  illuminated  by  our  instrument.  The  radiative  transport  equation  is  forced  to 
assume  a  perfect  point  source  for  illumination  -  for  measurements  over  long  distances  with  samples 
with  high  scattering  coefficients  and  fiber-launched  illumination  in  contact  with  the  sample,  the 
effective  angle  of  illumination  is  very  large  and,  after  travel,  many  scattering  events  will  occur, 
making  a  point-source  illumination  approximation  very  appropriate.  In  this  experiment,  a  collimated 
beam  is  used  as  the  source,  meaning  that  the  illumination  is  highly  directionalized  within  the  sample. 
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longer  focal  length 


The  difference  between  these  two  models 
is  illustrated  in  Figure  11a,  which  compares 
the  photon  fluence  at  four  points  on  a 
homogeneous  phantom  as  simulated  in 
Nirfast  and  measured  using  the 
experimental  setup.  It  is  immediately 
apparent  that  these  two  models  are  not 
equivalent.  With  a  point  source 
illumination,  photon  density  at  any  position 
in  the  phantom  is  directly  proportional  to 
the  distance  from  the  source.  Collection  at 
45°  relative  to  the  source  is  more  intense 
than  90°,  etc.  In  the  experimental  setup, 
fluence  is  greatest  at  a  collection  point  180° 
relative  to  the  source,  which  is  the  exact 
opposite  of  the  simulation.  In  this 
particular  experimental  setup,  the  source 
directionality,  scattering  level,  and 
measurement  distance  are  such  that  the 
assumption  of  photon  diffusion  behavior 
does  not  hold. 


d) 


Figure  12)  a)  Schematic  indicating  a  longer  detector  focal 
length  and,  in  turn,  a  larger  collection  area  by  the  detector, 
b)  PTFE  to  Intralipid  Raman  signal  ratio  averaged  over 
all  fibers.  Due  to  the  larger  area  of  collection,  the  peak 
corresponding  to  the  target  passing  the  detector  is  broader 
than  the  peak  corresponding  to  the  target  passing  the 
source,  c)  Fiber  optic  center  of  gravity  is  much  more 
varied  when  the  area  of  collection  is  increased,  possibly 
leading  to  more  angular  information  during  collection,  d) 
Nirfast  reconstruction  using  10  individual  fibers  spanning 
the  length  of  the  fiber  pattern  per  rotation  angle.  One 
target  is  predicted  with  no  shadowing. 

another  in  both  models.  As  such,  it  can  be  argued  that, 
around  180°  relative  to  the  source,  a  point  source 
experimental  setup's  collimated  beam  illumination. 


This  presents  an  interesting  problem  when 
working  with  measurements  where  the 
source  and  detectors  are  positioned  at 
multiple  sets  of  angles.  The  structure  of 
illumination  within  the  specimen  is 
fundamentally  different  between 
simulation  and  experiment  to  the  point 
where  the  radiative  transport  equation  fails 
to  describe  the  system  when  a  wide  range 
of  source-detector  angles  are  utilized. 
Figure  lib  shows  the  illumination 
structure  for  a  point  source  as  determined 
by  the  radiative  transport  equation  and 
Figure  11c  shows  an  example  of  a 
collimated  laser  source  as  determined  by  a 
ray  tracing  simulation.  While  the  fluence  in 
regions  between  the  source  and  detectors 
at  angles  smaller  than  180°  is  highly 
dependent  on  the  choice  of  model,  the 
fluence  at  detector  positions  with  angles 
very  close  to  180°  are  very  similar  to  one 
for  detector  positions  close  to  and  centered 
approximation  may  be  suitable  for  the 
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In  order  to  collect  at  a  number  of  angles  near  to  180°  from  the  source,  the  fiber  pattern  was  focused 
with  a  lens  of  length  200  mm,  significantly  increasing  the  size  of  the  fiber  pattern  as  projected  onto 
the  phantom  (as  illustrated  in  Figure  12a)  to  nearly  4  mm  across.  In  this  manner,  the  angles  formed 
between  the  source  and  detectors  ranged  from  0-10.4  degrees.  Upon  measuring  the  phantom  with 
this  setup  and  summing  the  responses  of  all  48  active  fibers,  the  resulting  response  shown  in  Figure 
12b  is  significantly  different  from  that  in  Figure  9a.  While  two  peaks  separated  by  approximately 
180°  are  observed,  one  peak  (corresponding  to  the  target  passing  the  detector  array)  is  significantly 
shorter  and  wider,  resulting  from  the  increased  spacing  of  the  detectors.  This  reasoning  is  confirmed 
in  Figure  12c,  in  which  the  standard  deviation  of  the  peak  position  for  all  fibers  is  increased  to  23.9. 
The  wider  spacing  of  the  fibers  implies  that  more  angular  information  is  collected.  Upon  performing 
a  Nirfast  reconstruction  using  10  individual  fibers  spaced  across  the  width  of  the  fiber  pattern  (and, 
as  a  result,  tenfold  the  number  of  data  points  as  compared  to  a  single  detector  at  180°  from  the 
source),  the  image  in  Figure  12d  is  generated.  A  single  target  is  predicted  with  minimal  'shadowing' 
artifacts  seen  from  collecting  at  angles  smaller  than  180°  from  the  source,  indicating  that  errors 
resulting  from  an  incorrectly  predicted  illumination  structure  are  reduced  when  source-detector 
angles  are  large  and  span  a  small  range. 

While  it  has  been  demonstrated  that  the  correct  number  of  targets  (and  an  approximation  of  target 
position)  can  be  recovered  using  this  large  fiber  pattern  setup,  the  validity  of  this  technique  cannot 
be  determined  unless  a)  the  separation  of  two  closely-spaced  targets  and  b)  the  accurate 
determination  of  the  positions  of  these  targets  can  be  demonstrated.  In  order  to  do  this,  three 
changes  were  made  to  the  experiment  summarized  in  Figures  8c  and  12a:  two  1/16"  diameter 
spheres  (instead  of  one  1/8"  diameter  sphere)  were  embedded  in  the  phantom  with  close  spacing 
(<=  sphere  diameter),  'fluorescence'  Nirfast  reconstructions  were  used  (as  to  separate  PTFE  and 
Intralipid  responses  and  normalize  each  to  a  simulated  blank),  and  micro-CT  was  performed  on  each 
phantom  in  order  to  determine  the  exact  target  positions. 


Figure  13)  a)  Micro-CT  reconstruction  of  the  phantom  with  two  PTFE  spheres,  b)  Nirfast  reconstruction,  c) 

Micro-CT  and  Nirfast  reconstruction  overlay. 

Figure  13a  shows  the  micro-CT  reconstruction  of  a  phantom  with  two  small  embedded  PTFE  spheres, 
illustrating  their  small  separation  relative  to  target  size.  Figure  13b  shows  the  Nirfast  reconstruction 
using  data  from  the  Raman  tomography  experiment.  Two  closely-spaced  targets  are  predicted  with 
some  shadowing  seen  in  the  center  of  the  reconstruction.  In  order  to  determine  how  accurate  the 
positions  and  spacing  of  these  predicted  targets  are,  Figure  13c  shows  an  overlay  of  the  micro-CT 
reconstruction  atop  the  Nirfast  reconstruction.  We  see  that  the  error  in  the  position  of  each  target  is 
less  than  1mm. 
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Figure  14)  a)  Micro-CT  reconstruction  of  phantom  with  two  targets  and  an  increased  (0.5%)  Intralipid 
content,  b)  Nirfast  reconstruction,  c)  Micro-CT  and  Nirfast  reconstruction  overlay. 

The  0.1%  Intralipid  phantoms  are  opaque  to  the  eye,  but  have  only  about  1/5  the  scattering  power  of 
some  epithelial  tissues.9  The  next  step  towards  translating  this  method  to  tissue  measurements  is  to 
increase  the  scattering  level  of  the  phantoms  fivefold  to  0.5%  Intralipid.  The  same  experiment  as 
shown  in  Figure  13  was  performed  on  a  phantom  consisting  of  0.5%  Intralipid.  Figure  14a  shows  the 
micro-CT  reconstruction  of  the  phantom,  Figure  14b  shows  the  Nirfast  reconstruction,  and  Figure 
14c  shows  the  overlay  of  the  two.  It  is  immediately  noticeable  that  the  positions  of  the  targets  are 
not  correct.  Although  the  angular  spacing  between  the  targets  is  correctly  reconstructed,  their 
positions  are  'projected'  onto  the  edge  of  the  phantom.  While  the  width  of  these  projections  are 
similar  to  the  widths  of  the  PTFE  spheres,  the  'depths'  of  these  targets  cannot  be  resolved  at  this  level 
of  scattering  using  the  current  instrumentation  and  reconstruction  methods. 

B:  Monte  Carlo  simulations 


Thus  far  it  has  been  shown  that  Raman  tomography  reconstructions  with  Nirfast  can  provide 
reasonable  size  and  position  information  for  buried  spheres.  As  experimental  parameters  become 
more  'extreme,'  such  as  a  fivefold  increase  in  the  scattering  power  of  the  medium,  reconstruction 
accuracy  is  rapidly  degraded.  Consideration  on  the  improvement  of  reconstruction  accuracy  begins 
with  an  exploration  of  factors  which  limit  reconstruction  quality.  Three  general  causes  for  poor 
reconstructions  that  will  be  discussed  here  include  not  acquiring  enough  total  data  (absolute 
undersampling),  not  acquiring  enough  unique  data  (spatial  undersampling),  and  misinterpreting 
acquired  data  (incorrect  modeling  of  system  response). 

Iterative  image  reconstruction  of  this  sort  often  falls  under  the  regime  of  ill-posed  inverse  problems, 
meaning  that  the  number  of  unknown  variables  to  be  solved  (the  number  of  pixels  in  the  resulting 
mesh)  is  larger  than  the  number  of  known  variables  (data  points  collected).  When  the  number  of 
measurements  performed  meets  or  exceeds  the  number  of  pixels  in  the  resulting  mesh,  a  direct 
solution  is  theoretically  possible.  In  any  case,  increasing  the  number  of  measurements  increases  the 
maximum  potential  reconstruction  quality.  At  the  same  time,  care  must  be  taken  to  ensure  that  each 
of  these  measurements  contributes  an  amount  of  unique  or  orthogonal  information  to  the  resulting 
data  set.  An  increase  in  the  absolute  number  of  measurements  cannot  improve  reconstruction 
quality  if  these  new  data  points  do  not  add  new  information  to  the  resulting  data  set.  With  the  end 
goal  of  determining  accurate  spatial  distributions,  the  combination  of  source  and  detector  sizes  and 
positions  needs  to  be  intelligently  chosen  such  that  no  region  of  the  specimen  is  poorly  probed  or 
under-sampled  as  compared  to  the  specimen  as  a  whole. 
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These  first  two  issues,  both  related  to  under-sampling,  can  be  addressed  to  a  certain  degree  without 
the  use  of  computer  modeling  or  an  exact  knowledge  of  the  system  response.  The  total  number  of 
measurements  to  acquire  can  be  chosen  based  on  the  desired  number  of  pixels  in  the  reconstruction 
mesh  (or  vice-versa),  whereas  choosing  a  wide  range  of  different  source-detector  configurations 
lends  itself  to  probing  a  number  of  different  regions  of  the  specimen.  Neither  of  these  approaches  is 
sufficient.  Without  a  strong  knowledge  of  both  which  regions  of  a  specimen  are  probed  in  a  given 
source-detector  configuration  and  the  collection  efficiency  of  each  configuration,  it  is  not  possible  to 
determine  the  effectiveness  of  a  given  instrument  setup.  With  such  knowledge,  an  ideal  instrument 
configuration  can  be  determined  and  the  maximum  reconstruction  quality  can  be  predicted. 

Knowledge  of  source-detector  responses  with  regards  to  a  given  specimen  (system-analyte  response) 
can  be  determined  through  the  computer  modeling  of  photon  fluence.  This  is  precisely  what  happens 
during  a  reconstruction  using  Nirfast:  photon  fluence  from  each  source  and  the  spatial  collection 
efficiency  of  each  detector  is  determined  through  radiative  transport  calculations  in  order  to 
quantify  the  sampled  region  of  each  source-detector  pairing.  As  mentioned  previously,  the  use  of 
radiative  transport  equations  in  the  diffusion  regime  requires  several  assumptions  about  both  the 
instrument  and  the  specimen.  One  such  assumption  is  that  all  sources  are  perfect  point  sources  and 
all  detectors  have  a  numerical  aperture  of  one  (equal  collection  efficiency  over  all  angles,  the  detector 
equivalent  of  a  point  source).  Additionally,  specimens  are  assumed  to  be  heavily  isotropically  light 
scattering  to  the  point  that  the  movement  of  photons  is  governed  by  the  diffusion  regime. 

While  these  assumptions  sufficiently  describe  DOT  instrumentation,  instruments  that  are  simply 
similar  to  DOT  will  deviate  from  these  assumptions  by  some  amount,  and  this  is  exemplified  by  the 
instrument  described  in  this  report.  The  illumination  is  not  a  point  source  directly  in  contact  with 
the  specimen,  but  is  instead  a  collimated  beam.  Detectors  are  also  not  in  contact  with  the  specimen 
such  that  each  detector's  response  changes  as  a  function  of  the  angle  it  makes  with  the  specimen 
surface.  Because  of  the  strong  directionalization  of  the  incoming  light,  the  scattering  power  of  the 
analyte,  and  the  small  distances  over  which  measurements  are  taken,  the  diffusion  approximation 
does  not  hold  (as  shown  in  Figure  11). 

Significant  deviations  from  the  idealized  diffusion  regime  have  two  main  impacts  on  Raman 
tomography  imaging  and  instrument  design:  optimal  instrumentation  configuration  and 
measurement  parameters  cannot  be  determined  if  accurate  specimen  illumination  and  light 
collection  cannot  be  modeled,  and  correct  reconstructions  cannot  be  obtained  if  reconstruction 
algorithms  are  misinterpreting  the  experimental  data  due  to  the  inability  to  accurately  determine  the 
specimen  region  probed  for  each  source-detector  pairing.  Correcting  for  these  differences  requires 
abandoning  the  diffusion  regime  and  its  required  assumptions  in  favor  of  a  modeling  technique  that 
allows  for  more  accurate  descriptions  of  instrument  sources,  detectors,  and  light-specimen 
interactions. 

Ray  tracing  through  an  iterative  Monte  Carlo  script  represents  a  solution  for  light  modeling  that 
accommodates  a  full  description  of  an  instrument's  sources,  detectors,  and  specimens  without 
operating  under  the  previously  discussed  assumptions  of  diffusion-regime  radiative  transport 
calculations.  Such  a  simulation  models  the  travel  of  individual  photons  as  they  are  generated  at  a 
source,  enter  a  specimen,  travel  through  the  specimen,  and  eventually  exit  the  specimen  or  become 
absorbed.  The  path  traveled  is  a  'random  walk'  governed  by  statistical  distributions  determined  by 
user-input  properties  of  the  instrument  and  specimen  (source  placement,  width,  and  angular 
distribution;  specimen  refractive  index,  mean  scattering  length,  and  mean  absorbance  length; 
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detector  placement,  width,  and  numerical  aperture).  By  simulating  and  recording  the  paths  of 
thousands  of  such  photons,  a  detailed  description  of  how  the  specimen  is  illuminated  and  the 
distribution  of  collected  photons  is  generated.  The  power  of  this  approach  is  that  it  is  relatively 
simple  and  straightforward  to  ray  trace  individual  photons  and  predict  their  behaviors  during 
interaction  with  the  specimen  via  Beer's  Law,  Snell's  Law,  Fresnel  Principle,  etc.  Modeling  thousands 
of  these  photons  does  not  require  additional  mathematics,  only  additional  computing  power.  As 
there  is  no  analytical  solution  for  the  radiative  transport  equation  outside  of  the  diffusion  regime, 
such  a  photon  modeling  approach  represents  a  reasonable  solution  to  more  accurate  modeling  of  a 
range  of  instruments  and  specimens. 

Monte  Carlo  approaches  are  not  without  drawbacks.  Modeling  a  single  photon  is  a  mildly 
computationally  intensive  iterative  process;  modeling  thousands  of  photons  for  each  source  and 
detector  can  increase  this  computational  burden  by  several  orders  of  magnitude.  While  there  have 
been  solutions  in  the  literature  to  increase  computational  efficiency,  including  modeling  'big'  photons 
which  undergo  multiple  absorbance  eventsi13]  and  the  use  of  graphics  processing  units  (GPUs)  for 
code  speedup, t14l  at  the  present  time  Monte  Carlo  methods  take  strictly  more  computing  power  than 
radiative  transport  calculations.  Additionally,  simulations  driven  by  underlying  random  processes 
do  not  converge  on  an  exact  solution.  Just  as  the  measurement  of  a  signal  containing  white  noise 
approaches  the  true  value  as  acquisition  time  is  increased,  Monte  Carlo  simulations  approach  an 
ideal  solution  as  more  photons  are  modeled,  but  there  will  always  be  variability  in  the  result. 
Reducing  this  variability  to  an  acceptable  level  through  increased  runtimes  is  another  consideration 
when  employing  these  methods.  If  these  computational  difficulties  can  be  handled,  these  methods 
can  provide  a  straightforward  and  fully  customizable  means  of  modeling  specimen  illumination  and 
performing  reconstructions  from  experimental  data. 


Feature 

Radiative  Transport 
+  Diffusion 

Monte  Carlo 

Photon  migration 

Diffusion  behavior 

No  diffusion 
assumption 

Source 

Point  sources 

Directionalized 
(laser)  illumination 

Detectors 

NAof  1 

NA<  1 

Table  1.  Comparison  of  the  three  main  differences  between  radiative  transport  (Nirfast)  and  Monte  Carlo 
simulations.  The  lack  of  an  analytical  solution  to  radiative  transport  equation  outside  of  the  diffusion  regime 
requires  that  photons  migration  mimics  diffusion  and  that  all  point  sources  and  detectors  are  idealized. 

Monte  Carlo  simulations  do  not  require  these  assumptions,  as  photons  can  migrate  outside  of  the  diffusion 
regime  and  sources  and  detectors  can  be  modeled  to  have  small  numerical  apertures. 

To  demonstrate  some  of  these  advantages,  a  Monte  Carlo  photon  ray  tracing  code  was  written  in- 
house  with  Matlab.  The  cylindrical  specimen  was  represented  by  a  circular  mesh  (similar  to  that 
used  in  the  Nirfast  reconstructions)  comprised  of  709  square  pixels  (as  compared  to  ~1750 
triangular  elements  in  the  Nirfast  reconstructions),  each  with  a  defined  mean  scattering  length,  mean 
absorbance  length,  refractive  index,  etc.  For  the  following  discussion,  fluence  modeling  was 
performed  on  a  simulated  homogeneous  circular  tissue  phantom  with  16  equally-spaced  sources  and 
detectors  spread  equally  over  the  surface.  This  is  often  referred  to  as  a  standard  'fan'  geometry  in 
DOT  literature!1]  For  the  Nirfast  simulations,  all  sources  were  treated  as  point  sources,  all  detectors 
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had  an  effective  numerical  aperture  of  one,  and  photon  migration  occurred  under  the  assumption  of 
the  diffusion  regime.  These  three  characteristics  are  required  for  Nirfast  simulations  and  diffusion- 
regime  radiative  transport  calculations  in  general.  In  contrast,  the  Monte  Carlo  simulations  utilized 
directionalized  laser  sources,  small  numerical  aperture  detectors,  and  a  Poissonian  distribution  to 
model  scattering  such  that  diffusion  behavior  was  not  assumed.  While  the  assumptions  required  by 
radiative  transport  calculations  have  worked  well  in  the  DOT  field,  the  small  specimen  size  and 
directionalized  laser  illumination  utilized  in  these  Raman  tomography  experiments  are  not  well 
described  by  these  assumptions.  These  differences  are  summarized  in  Table  1.  The  images  from  the 
Monte  Carlo  simulation  shown  below  look  significantly  more  'pixelated'  than  their  Nirfast 
counterparts  because  the  data  are  not  heavily  interpolated  before  being  presented  to  the  user. 


Figure  14.  Simulated  sampling  regions  of  a  homogeneous  tissue  phantom  for  two  detectors  separated  by 
157.5  degrees,  a)  Nirfast  simulation  with  a  point  source,  idealized  detector,  and  diffusion-regime  photon 
migration,  b)  Monte  Carlo  simulation  modeling  a  directionalized  laser  source,  a  small  numerical  aperture 

detector,  and  non-diffusion  photon  migration. 

Visualizing  the  sampling  region  of  a  single  source-detector  pair  using  the  adjoint  source 
approximation  illustrates  the  differences  between  point  and  directionalized  sources  and  detectors. 
The  adjoint  source  approximation,  employed  by  both  Nirfast  and  the  Monte  Carlo  simulation,  states 
that  if  a  source  and  a  detector  switch  positions  the  detected  power  will  remain  the  same  as  the 
detected  photons  must  travel  through  the  same  region.  Simulating  the  fluence  from  a  source  and  the 
fluence  from  a  detector  (as  though  it  were  a  source)  and  multiplying  these  two  regions  together 
shows  the  mutual  'sampling'  region  of  the  two.  This  approximation  holds  true  as  long  as  scattering  is 
assumed  to  be  isotropic  or  if  the  optical  process  bridging  the  two  regions  (spontaneous  Raman 
scattering)  is  isotropic.  Figure  14  depicts  a  homogeneous  circular  mesh  with  a  source-detector  pair 
at  an  angle  of  157.5°.  Using  the  adjoint  source  approximation,  the  effective  sampling  region  for  both 
a  point  source  with  a  'point'  detector  and  a  directionalized  laser  source  with  a  less-than-unity 
numerical  aperture  detector  are  shown.  The  Monte  Carlo  simulation  shows  that  the  majority  of  the 
collected  signal  originates  from  the  bulk,  while  in  the  diffusion  regime  the  most  intense  signal 
originates  from  the  phantom's  surface.  Quantitatively,  the  sampling  power  distribution  of  these  two 
regions  differs  by  approximately  23%.  These  differences  in  sampling  region  are  more  pronounced 
when  all  source-detector  combinations  are  taken  into  account,  representing  the  sampling  region  for 
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Figure  15)  Simulated  sampling  regions  from  all  240  source-detector  pairs  for  a  homogeneous  circular  tissue 
phantom.  A)  In  the  diffusion  regime  (very  high  scattering,  large  sampling  distances),  the  majority  of  the 
collected  signal  originates  from  near  the  phantom’s  surface.  B)  In  a  lower-scattering  regime  with 
directionalized  illumination,  there  is  a  strong  preference  for  signal  originating  from  the  bulk.  An  X- shaped 

artifact  is  present  due  to  the  coarse  pixel  mesh. 

an  entire  experiment.  The  results  from  summing  the  sampling  regions  for  all  240  source -detector 
combinations  are  shown  in  Figure  15.  For  the  Monte  Carlo  simulation  representing  our  Raman 
tomography  setup,  the  majority  of  the  collected  signal  originates  from  the  phantom's  bulk,  while  a 
diffusion  regime  behavior  results  in  a  collection  preference  for  the  phantom's  surface.  A  strong 
preference  for  bulk  signal  agrees  with  previous  Monte  Carlo  simulations  regarding  to  transmission 
Raman  measurements!11] 


Figure  16)  Monte  Carlo  simulations  can  be  used  for  the  fluence-modeling  ‘forward’  step  in  an  attenuation- 
based  reconstruction,  a)  A  homogeneous  mesh  with  an  anomalous  region  of  absorbance  that  is  lOx  higher 
than  the  bulk,  b)  After  five  iterations,  a  Monte  Carlo-based  reconstruction  algorithm  can  identify  the  region 

of  the  anomaly,  albeit  with  artifacts. 

Fluence  modeling  simulations  through  Monte  Carlo  methods  can  also  be  used  for  reconstruction 
purposes.  Using  the  same  mesh,  a  region  of  anomalous  absorbance  was  added  (Figure  16),  simulated 
experimental  data  was  generated  by  Monte  Carlo  fluence  modeling,  and  the  data  was  supplied  to  a 
reconstruction  algorithm  utilizing  the  same  Monte  Carlo  method  and  homogeneous  absorbance 
properties  for  a  starting  guess.  It  should  be  noted  that  this  reconstruction  type  is  fundamentally 
different  from  those  discussed  in  the  rest  of  this  report,  as  Raman/fluorescence  reconstructions 
investigate  a  photon  generation  process  and  require  measured  values  for  both  the  Raman  signal  and 
the  fundamental  frequency,  while  an  'absorbance'  reconstruction  is  attenuation  based  and  only 
utilizes  one  data  point  per  source-detector  pair.  Even  so,  this  demonstrates  that  Monte  Carlo  fluence 
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simulations  have  sufficient  repeatability  to  be  useful  in  reconstruction  algorithms.  In  order  to 
maximize  this  utility,  processing  times  need  to  be  decreased  significantly  through  the  use  of  graphics 
processing  units  (GPUs). 

These  simulation  comparisons  suggest  that  while  diffuse  optical  tomography  and  Raman  tomography 
are  similar  in  form  and  function,  they  can  differ  drastically  in  analyte  response.  Radiative  transport 
calculations  have  been  an  invaluable  tool  in  the  DOT  field,  as  the  results  are  highly  reproducible  and 
accurately  model  those  experiments.  Upon  moving  to  lower  scattering  levels,  smaller  sample  sizes, 
and  the  use  of  laser  sources  and  small-NA  detectors,  the  diffusion-regime  radiative  transport 
equation  becomes  a  less  accurate  model.  In  order  to  develop  the  field  of  transmission  Raman 
tomography,  elegant  fluence  modeling  simulations  allowing  for  the  accurate  modeling  of  arbitrary 
sources,  detectors,  and  sample  geometries  will  need  to  be  realized.  Monte  Carlo  simulations  bolstered 
by  the  processing  speed  of  graphics  processing  units  represent  one  possible  means  to  this 
development. 

Work  aimed  at  Raman  tomographic  reconstruction  is  now  underway  using  the  Monte  Carlo  approach 
towards  aim  2.3  for  year  2  of  this  research. 

An  additional  point  worth  adding  to  this  report  is  in  our  collection  scheme.  Through  the  research 
conducted  in  the  above  experiments  and  through  modeling  we  determined  a  fiber  placed  in  contact 
with  the  tissue  surface  should  provide  the  highest  collection  efficiency.  The  maximum  collection 
efficiency  occurs  with  the  fiber  in  contact  with  the  breast  tissue.  The  collection  efficiency  of  a  fiber 
optic  probe  (the  absolute  number  of  photons  that  reach  the  CCD  camera)  increases  when  the  fiber  is 
in  direct  contact  with  the  tissue  as  opposed  to  collection  through  a  distance  of  air  via  a  focusing  lens. 
The  two  phenomena  that  govern  this  effect,  Fresnel  reflection  and  numerical  aperture,  are  primarily 
functions  of  the  refractive  indices  of  tissue  (n  ~=  1.3-1. 5)  t15l,  silica  glass  (n  ~=  1.5-1. 6),  and  air  (n  = 
1).  The  magnitude  of  Fresnel  reflections  at  the  boundary  of  two  materials  decreases  as  the  difference 
between  the  materials'  refractive  indices  (An)  approaches  zero.  First,  eliminating  a  'layer'  of  air 
between  tissue  and  the  fiber  probe  reduces  the  number  of  interfaces  at  which  reflections  occur. 
Furthermore,  the  difference  in  refractive  index  for  a  tissue-fiber  boundary  (An  ~=  0  to  0.3)  is  smaller 
than  either  a  tissue-air  boundary  (An  ~=  0.3-0. 5)  or  an  air-fiber  boundary  (An  ~=  0.6).  In  addition, 
the  effective  numerical  aperture  (or  collection  angle)  of  a  fiber  optic  or  lens  increases  as  the 
refractive  index  interfacing  with  the  specimen  increases.  Just  as  an  oil-immersion  microscope 
objective  has  a  larger  numerical  aperture  than  a  standard  objective  that  operates  in  airt16!, 
'immersing'  the  tissue  in  a  refractive  index  n~=1.5-1.6  medium  by  directly  interfacing  the  fiber  and 
the  tissue  increases  both  the  collection  angle  and  the  total  number  of  photons  coupled  into  the  fiber. 
Therefore  in  our  instrument's  design  we  have  included  breast  paddles  with  slots  for  positioning 
fibers  so  that  we  are  able  place  individual  fiberoptics  in  direct  contact  with  the  skin's  surface. 
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KEY  RESEARCH  ACCOMPLISHMENTS: 

•  We  have  become  familiar  with  clinical  practices  associated  with  breast  cancer  screening, 
breast  cancer  histology,  how  breast  tissue  is  sampled  and  biopsied  for  diagnosis,  and  breast 
cancer  pathology. 

•  We  have  built-up  a  comprehensive  database  consisting  of  micrographs  and  IR  spectral 
images  for  identifying  breast  tissue  histology  and  tissue  chemistry.  This  database  contains 
tissue  cases  diagnosed  as  hyperplasia,  dysplasia,  malignant  carcinomas,  and  normal  tissue. 

•  Using  this  database  we  have  identified  several  locations  on  each  biopsy  to  acquire  Raman 
measurements 

•  We  have  acquired  Raman  measurements  on  all  biopsy  sections  as  well  as  several  surgical 
resections. 

•  We  have  identified  Raman  spectral  bands  that  can  be  used  for  distinguishing  between 
different  tissue  types. 

•  We  have  used  the  information  obtained  from  the  point  Raman  measurements  to  achieve  cell- 
level  contrast  in  Raman  spectral  images 

•  We  have  used  clinical  observations  and  interaction  with  clinicians  to  develop  a  conceptual 
design  for  a  Raman  Tomography  instrument  aimed  at  breast  cancer  screening. 

•  We  have  put  together  autoCad  files  and  have  begun  construction  on  the  instrument 

•  We  have  evaluated  existing  diffuse  optical  tomography  algorithms  for  Raman  tomography 
and  have  adapted  a  Monte  Carlo  framework  for  Raman  tomographic  reconstruction 

•  We  have  investigated  through  modeling  and  experimentation  different  instrument 
configurations  and  have  implemented  the  transmission  fan  style  configuration  into  our 
instrument  design 

•  We  have  determined  that  a  fiber  in-contact  with  tissue  should  give  us  the  maximum 
collection  efficiency 

REPORTABLE  OUTCOMES: 

•  A  manuscript  detailing  high  resolution  IR  images  is  included  in  the  appendix  of  this  report. 
This  manuscript  is  indirectly  a  result  of  this  work  in  that  we  were  trying  to  obtain  high 
resolution  IR  images  for  identifying  the  tissue  chemistry  of  biopsies.  This  information  was 
employed  in  determining  regions  for  Raman  point  measurements.  This  manuscript  has  been 
accepted  for  publication  in  the  Journal  of  Applied  Spectroscopy. 

•  A  manuscript  detailing  'data  mining'  and  visualization  of  hyperspectral  data  is  included  in  the 
appendix  of  this  report.  This  manuscript  was  the  result  of  looking  for  spectral  correlations  in 
different  tissue  types  and  reports  an  approach  for  efficient  visualization  of  such  a  data  set. 

The  manuscript  is  currently  under  review  at  the  Journal  of  Transactions  on  Visualization  and 
Computer  Graphics. 

•  We  are  working  on  a  manuscript  reporting  cell  type  classification  through  Raman 
spectroscopy  which  will  detail  the  spectral  bands  correlated  with  the  different  cell  types. 
From  the  year  one  work  we  have  identified  Raman  bands  such  as  the  1045  cm1  Raman  shift, 
the  828cm1  O-P-O  stretch,  and  the  amide  III  bands  which  can  be  used  as  chemical  contrast 
for  distinguishing  different  cell  types. 

•  We  are  working  on  a  manuscript  reporting  a  comparison  in  diffuse  optical  tomography 
reconstruction  for  Raman  spectroscopy 

•  Oct.  2011-  Poster  Presentation  46th  Midwest  Regional  ACS  Meeting,  Saint  Louis  MO 

•  Feb.  2012-  Oral  Presentation  at  the  University  of  Michigan,  Ann  Arbor 

•  March  2012  -  Poster  presentation  at  the  University  of  Illinois  Chicago  Cancer  Center  Forum 

24 


CONCLUSION:  We  have  successfully  completed  the  work  described  in  the  approved  statement  of 
work  for  year  1  of  this  research  program.  The  results  from  this  first  year  have  forged  a  nice 
foundation  for  the  work  we  have  planned  in  year  2.  We  now  have  a  working  relationship  with  local 
clinicians  and  have  become  familiar  with  clinical  practices  associated  with  breast  cancer  screening, 
breast  cancer  histology,  how  breast  tissue  is  sampled  and  biopsied  for  diagnosis,  and  breast  cancer 
pathology.  We  have  built-up  a  comprehensive  database  consisting  of  micrographs  and  IR  spectral 
images  with  identified  breast  tissue  histology  and  tissue  chemistry.  This  database  has  tissue  showing 
cases  diagnosed  as  hyperplasia,  dysplasia,  malignant  carcinomas,  and  normal  tissue.  Using  this 
database  we  have  identify  several  locations  on  each  biopsy  to  acquire  Raman  measurements  and 
have  acquired  Raman  measurements  on  biopsy  sections  as  well  as  several  surgical  resections. 
Through  careful  analysis  we  have  identified  Raman  spectral  bands  that  can  be  used  for  distinguishing 
between  different  tissue  types  and  have  applied  those  Raman  spectral  bands  to  achieve  cell-level 
contrast  in  Raman  spectral  images.  Additionally,  we  have  used  clinical  observations  and  interaction 
with  clinicians  to  develop  a  conceptual  design  for  a  Raman  Tomography  instrument  aimed  at  breast 
cancer  screening.  From  this  conceptual  design  we  put  together  autoCad  files  and  have  begun 
construction  on  a  prototype  instrument  for  conducting  Raman  tomography  on  breast  tissue.  For 
Raman  tomographic  reconstruction,  we  have  evaluated  existing  diffuse  optical  tomography 
algorithms  for  specifically  for  Raman  measurements  and  have  adapted  a  Monte  Carlo  framework  for 
our  Raman  tomographic  reconstruction.  Through  modeling  and  experimentation  we  characterized 
different  instrument  configurations  and  have  implemented  the  transmission  fan  style  configuration 
into  our  prototype  instrument  design  which  uses  fiber-optics  in-contact  with  tissue  to  give  the 
maximum  collection  efficiency.  With  these  accomplishments  we  are  closer  to  achieving  prototype 
instrument  capable  of  Raman  tomography  in  breast  tissue.  We  are  excited  to  continue  this  work  in 
year  2. 
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Designing  Transfer  Functions  for  Interactive 
Data  Mining  of  Hyperspectral  Images 

David  Mayerich,  Member,  IEEE,  Michael  Walsh,  Matthew  Schulmerich,  Rohit  Bhargava 


Abstract — Spectroscopy  is  used  in  several  fields  to  acquire  chemical  information  about  a  substance  for  identification  or  analysis. 
Recent  advances  allow  the  acquisition  of  hyperspectral  imagery  for  biomedical  data.  This  has  promising  implications  in  the  fields 
of  biotechnology  and  medicine,  where  quantitative  analysis  of  tissue  can  be  performed  by  directly  measuring  spatially-resolved 
chemical  information.  At  present,  however,  these  data  sets  require  a  significant  amount  of  human  intervention  to  identify  spectral 
features  for  use  in  classification.  We  propose  a  method  for  designing  transfer  functions  for  hyperspectral  images.  This  method 
allows  researchers  to  interactively  adjust  parameters  used  to  systematically  change  the  input  spectra  in  order  to  find  metrics  that 
can  be  used  to  classify  features  in  the  images. 

Index  Terms — microscopy,  spectroscopy,  transfer  functions,  rendering. 

- ♦ - 


1  Introduction 

Hyperspectral  data  is  composed  of  a  series  of  mea¬ 
surements,  taken  across  many  wavelengths  in  the 
electromagnetic  spectrum.  Typically,  materials  exhibit 
a  characteristic  spectral  signature,  that  is  specifically 
indicative  of  their  molecular  composition.  While  the 
optical  spectrum  is  commonly  used  for  molecular 
studies,  other  forms  of  analysis  of  a  wide  //spectrum// 
of  properties  has  led  to  the  development  of  a  class 
of  techniques  for  molecular  analysis.  Some  of  the 
most  common  examples  include  (a)  vibrational  spec¬ 
troscopy,  which  is  used  to  identify  chemical  species 
based  on  the  electromagnetic  recording  of  oscillations 
of  molecular  bonds  present  in  the  sample,  (b)  mass 
spectrometry,  which  is  used  to  identify  compounds 
based  on  the  mass-to-charge  ratio  of  their  constituent 
particles,  and  (c)  nuclear  magnetic  resonance  spec¬ 
troscopy,  which  allows  the  characterization  of  nuclei 
based  on  their  behavior  within  a  magnetic  field.  These 
techniques  are  often  used  in  the  fields  of  biomedicine 
and  forensics  to  determine  the  identity  of  chemi¬ 
cal  compounds.  While  each  has  its  advantages  and 
specific  areas  of  applicability,  optical  techniques  are 
most  useful  for  label-free  microscopy  with  micron- 
scale  resolution  and  significant  molecular  detail.  Here, 
we  use  data  from  these  techniques  as  a  demonstrative 
example  for  our  developed  methods. 

Recent  advances  in  optical  detector  technology  al¬ 
low  the  rapid  acquisition  of  spectral  signals  that  are 
spatially  specific,  producing  multidimensional  images 
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that  can  be  extensive  in  size  (several  tens  of  gi¬ 
gabytes).  These  techniques  have  shown  promise  in 
biomedicine  for  cell  type  classification  [14]  and  cancer 
analysis  [21].  However,  the  interpretation  of  a  spectral 
signature  is  a  complex  task,  requiring  a  significant 
amount  of  pre-processing  before  identifying  spectral 
features  that  correspond  to  structural  and  chemical 
information. 

While  the  acquisition  of  data  is  relatively  rapid, 
there  are  limited  options  at  present  to  assist  the  prac¬ 
titioner  in  analyzing  data.  In  this  paper,  we  demon¬ 
strate  techniques  for  designing  transfer  functions  that 
allow  a  user  to  interactively  explore  hyperspectral 
images.  We  demonstrate  these  techniques  on  data  sets 
imaged  using  vibrational  spectroscopy  techniques, 
with  a  particular  focus  on  mid-infrared  spectroscopic 
imaging,  where  structural  and  chemical  information 
are  particularly  difficult  to  separate. 

2  Vibrational  Spectroscopy 

Vibrational  spectroscopy  is  used  to  identify  a  chemical 
compound  based  on  molecular  bond  vibration.  Atoms 
that  compose  a  molecular  bond,  for  example  the  one 
between  carbon  and  hydrogen  (C-H),  are  not  static 
but  oscillate  at  specific  frequencies.  Any  chemical 
compound  and  its  molecular  organization  can  be 
identified  by  determining  the  type  and  number  of 
resonant  bonds.  This  information  can  be  determined 
by  examining  the  interaction  of  electromagnetic  ra¬ 
diation  with  the  specimen.  Two  common  methods  for 
molecular  spectroscopy  are  mid-infrared  spectroscopy 
and  Raman  spectroscopy  (Figure  1). 

2.1  Mid-Infrared  Spectroscopy 

Mid-infrared  spectroscopy  is  used  to  identify  molec¬ 
ular  bonds  by  determining  the  fraction  of  incident 
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(a)  Mid-Infrared  Absorbance  Spectroscopy 


(b)  Raman  Spectroscopy 


Fig.  1.  Mid-infrared  and  Raman  spectroscopy,  (a)  Mid- 
infrared  spectroscopy  measures  the  light  absorbed 
due  to  resonance  with  specific  vibrational  modes  in  a 
molecule.  When  transmitted  through  a  material  con¬ 
taining  a  CH2  functional  group,  as  shown  in  this  exam¬ 
ple,  the  output  signal  is  lower  intensity  than  the  incident 
light  when  ^  =  2853cm-1.  (b)  Raman  spectroscopy 
measures  the  same  vibrational  modes  by  applying 
an  input  wavelength  vim  A  very  small  fraction  of  the 
incident  light  is  inelastically  scattered,  resulting  in  a 
shift  from  the  incident  wavelength  by  ±2853cm-1  for 
CH2. 

light  absorbed.  Broadband  mid-infrared  electromag¬ 
netic  radiation  is  transmitted  through  the  sample  and 
the  absorbed  radiation  at  each  wavelength  is  quan¬ 
tified  (Fig.  la).  The  most  common  instrumentation 
for  measuring  absorption  spectra  is  Fourier  Transform 
Infrared  (FT-IR)  Spectroscopy  [16].  In  this  instrumen¬ 
tation,  and  uses  an  interferometer  is  employed  to 
encode  optical  frequencies  in  time  and  the  recorded 
data  are  decoded  using  a  a  Fourier  transform.  The 
resulting  data  consist  of  the  light  transmitted  as  a 
function  of  frequency.  By  measuring  the  spectrum 
with  and  without  a  sample,  hence,  the  absorbance 
of  the  sample  can  be  calculated.  Molecular  bonds  are 
therefore  identified  by  finding  peaks  in  the  resulting 
absorbance  spectrum  with  the  frequency  domain  be¬ 
ing  generally  recorded  in  units  of  wavenumber  (cm-1). 
This  technique  has  been  used  extensively  in  forensics 
[3],  chemistry,  and  biology. 

Recent  advances  in  FT-IR  imaging  instrumentation 
allow  the  use  of  focal  plane  arrays  (FPAs)  for  acquir¬ 
ing  spatially-resolved  mid-infrared  absorbance  spec¬ 
tra  at  high  speeds  [23].  This  allows  the  facile  collection 
of  hyperspectral  images,  where  each  pixel  provides 
the  corresponding  spatially  resolved  absorption  as  a 
function  of  wavenumber. 

2.2  Raman  Spectroscopy 

Raman  spectroscopy  provides  a  method  for  measur¬ 
ing  the  vibration  of  molecular  bonds  using  monochro¬ 
matic  light  [24].  Conventional  Raman  spectroscopy 
has  applications  in  fields  including  pharmaceuticals 
[12],  biomedicine  [6],  art  restoration  [1],  and  forensics 


[28].  Raman  spectroscopy  is  conventionally  performed 
by  focusing  narrow-band  (laser)  light  onto  a  sample. 
The  vast  majority  of  photons  are  elastically  scattered, 
maintaining  the  frequency  of  the  incident  light.  A 
small  number  of  photons  (1  in  10  million)  loose  energy 
equal  to  the  frequency  of  the  molecular  vibrations 
present  in  the  sample  and  therefore  are  a  way  of 
probing  the  molecules  that  are  present  (Fig.  1).  Raman 
images  are  conventionally  acquired  using  a  micro¬ 
scope  in  point  measurement  mode.  A  notch  filter  is 
used  to  reject  the  incident  wavelength  while  scattered 
light  is  dispersed  onto  a  CCD  array  using  a  grating. 


Traditional  methods  for  determining  the  chemical 
composition  of  a  specimen  rely  on  collecting  a  single 
spectrum  from  a  homogeneous  sample.  This  tech¬ 
nique,  known  as  point  spectroscopy,  was  the  standard 
in  vibrational  spectroscopy  due  to  both  the  cost  of 
infrared  (IR)  detector  arrays  and  the  large  time  needed 
to  acquire  Raman  spectra.  The  increasing  availability 
of  imaging  systems  now  permit  the  analyses  of  non- 
homogeneous  samples,  where  structural  changes  ac¬ 
company  changes  in  chemical  composition.  Chemical 
density  and  sample  thickness  play  a  role  in  light  ab¬ 
sorption.  These  structural  changes  also  result  in  elas¬ 
tic  scattering,  which  induces  wavelength-dependent 
changes  in  the  measured  spectra  [5],  [9],  [10].  In 
addition,  the  extraction  of  chemical  information  from 
an  absorption  spectrum  is  a  complex  task,  often  re¬ 
quiring  prior  knowledge  about  the  specimen  type  and 
possible  composition  [11].  In  this  section,  we  briefly 
describe  how  sample  structure  affects  a  spectrum 
as  well  as  techniques  for  de-coupling  structural  and 
chemical  information. 


The  refractive  index  of  a  material  can  be  expressed  as 
the  complex  function 

n{v)  =  r](u)  +  in{v)  (1) 

where  the  index  of  refraction  n  is  the  sum  of  the 
real  component  r\,  referred  to  as  the  phase  speed,  and 
the  imaginary  component  k,  referred  to  as  the  ex¬ 
tinction  coefficient.  Both  of  these  components  change 
as  a  function  of  wavelength.  The  real  component  rj 
helps  to  describe  the  change  in  light  direction  as 
it  crosses  boundaries  of  different  materials,  and  is 
a  key  component  in  ray-tracing  algorithms  that  in¬ 
corporate  light  refraction  [15],  [37].  Light  absorption 
due  to  the  resonant  vibration  of  molecular  bonds  is 
integrated  into  n.  However,  both  of  these  expressions 
provide  important  information  about  the  structure 
and  composition  of  a  specimen.  The  measurements 
performed  using  mid-infrared  spectroscopic  imaging 
provide  indirect  measurements  of  both  and  rj{v) 
through  the  measured  absorbance  A{v). 


3  Data  Mining 


3.1  Refractive  Index 
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Fig.  2.  Scattering  of  incident  light  through  non- 
homogeneous  samples.  Two  prominent  contributions 
to  scattering  are  (a)  edge  effects  and  (b)  scattering 
due  to  morphologic  features  in  the  sample.  Scattering 
causes  spatially  varying  and  wavelength-dependent 
fluctuations  in  transmission  intensity. 

These  measurements  are  tightly  coupled  and  must 
be  separated  in  order  to  determine  useful  information 
about  the  sample.  However,  the  separation  of  these 
components  is  a  difficult  problem  that  requires  man¬ 
ual  exploration  of  large  hyperspectral  images  as  well 
as  some  prior  knowledge  about  the  content  of  these 
images.  The  researcher  must  first  separate  spectral 
features  caused  by  light  scattering  via  micro-scale 
interactions  with  cells,  membranes,  and  other  edges  in 
a  biological  tissue  sample.  These  features  can  then  be 
removed  from  the  spectrum,  leaving  only  components 
of  absorption. 

3.2  Elastic  Scattering 

Light  that  is  transmitted  through  non-homogeneous 
samples  is  subject  to  scattering  at  boundaries  between 
materials  exhibiting  different  phase  speeds  r\.  These 
effects  are  prominent  in  mid-infrared  spectroscopy, 
where  the  re-direction  of  light  out  of  the  collection 
optics  of  the  microscope  can  appear  as  absorbance.  In 
addition,  scattering  effects  vary  based  on  wavenum¬ 
ber  due  to  the  wavelength-dependent  phase  speed 
r\{v')  (Sec.  3.1). 

These  effects  are  particularly  prominent  at  edges 
and  through  structures  at  scales  near  the  wavelength 
of  the  incident  light  (Fig.  2).  The  latter  has  been  mod¬ 
eled  using  the  formalism  of  resonant  Mie  scattering, 
which  describes  light  scattering  through  a  sphere. 
This  form  of  scattering  appears  as  a  low-frequency 
oscillation  in  the  absorbance  spectrum.  This  oscillation 
is  visible  when  imaging  cells  and  cell  nuclei.  Methods 
have  been  proposed  for  adjusting  spectra  based  on 
Mie  theory  [2],  however  these  techniques  are  time 
consuming  and  would  not  provide  interactive  feed¬ 
back.  In  addition,  Mie  theory  is  not  generally  appli¬ 
cable  to  scattering  effects  and  the  prior  information 
required  for  these  estimations  is  not  always  available. 
Finally,  these  models  do  not  account  for  microscopy 
optics. 


3 

The  removal  of  scattering  effects  can  be  accom¬ 
plished  through  other  models  and  is  termed  baseline 
correction.  For  example,  a  second  model  is  more 
commonly  used  in  simply  subtracting  a  piecewise- 
linear  function.  We  have  found  that  this  first-order 
approximation  for  baseline  correction  is  effective  for 
visualization.  We  allow  the  user  to  specify  points 
where  chemical  absorption  is  assumed  to  be  zero.  We 
then  dynamically  adjust  a  histogram  of  the  spectra  to 
reflect  these  parameters.  This  allows  the  user  to  lever¬ 
age  prior  knowledge  about  the  chemical  composition 
of  the  specimen,  while  providing  interactive  feedback 
that  allows  them  to  refine  their  selection. 

3.3  Beer-Lambert  Law 

After  contributions  due  to  scattering  are  identified 
and  removed,  the  remaining  spectral  features  corre¬ 
spond  to  the  absorbance  A  of  the  specimen.  These 
features  identify  light  absorbance  due  to  resonant 
molecular  vibrations  in  the  sample  as  a  function  of 
wavelength.  However,  physical  characteristics  of  the 
sample  still  affect  the  value  of  A  in  the  form  of  tissue 
thickness  and  density.  In  the  case  of  transmitted  light, 
these  relations  are  expressed  using  the  Beer-Lambert 
law: 

T  =  e-K(N  (2) 

where  £  is  the  path  length  through  the  specimen 
(thickness),  N  is  the  density  of  a  chemical  compound, 
and  hi  is  the  extinction  coefficient.  By  applying  this 
technique  to  absorbance  spectra  and  specifying  the 
extinction  coefficient  as  a  function  of  wavelength,  the 
Beer-Lambert  law  can  be  represented  as: 

A{p)  =  k(D)£N  (3) 

Note  that  there  are  no  methods  for  separating  the 
path  length  £  and  density  N  in  the  general  case. 
However,  in  many  cases  the  function  n{p)  can  be 
estimated  using  a  reference  band  vr  within  the  spec¬ 
trum.  One  example  of  a  useful  reference  is  the  Amide 
I  absorption  peak  often  found  in  biological  tissue 
samples  at  vr  ^1650cm-1.  By  normalizing  Amide 
I  absorption  to  be  uniform  across  chemical  species 
(hv(Pr)  =  1.0),  the  combined  contributions  of  path 
length  and  density  can  be  estimated  in  terms  of  the 
relative  protein  composition  using  A{pr)  —  £N  and 
the  extinction  coefficient  as: 


A  simple  example  of  using  a  reference  wavelength 
to  separate  chemical  components  is  shown  in  Figure 
3.  While  a  single  reference  feature  is  useful,  it  is  not 
generally  applicable  to  most  samples,  therefore  sev¬ 
eral  references  are  often  used  for  classifying  multiple 
chemical  species. 
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(a)  (b) 


Fig.  3.  Structure  and  chemical  features  of  an  ab¬ 
sorbance  spectrum,  (a)  A  specimen  is  represented 
as  a  collection  of  two  tissue  types  of  varying  den¬ 
sity/thickness  placed  on  a  barium  fluoride  substrate 
(top)  with  the  corresponding  density  (bottom),  (b)  The 
IR  spectra  from  regions  C0,  C±,  and  C2  are  domi¬ 
nated  by  stattering.  (c)  The  spectra  are  shown  after 
specifying  baseline  points  (arrows).  Thickness  can  be 
determined  at  the  reference  peak,  (d)  Normalizing  the 
spectra  using  a  reference  allows  separation  of  the 
chemical  compounds  with  a  single  feature. 


4  Transfer  Function  Design 

Transfer  functions  are  commonly  used  in  volume 
visualization  [13],  [32]  to  assign  color  and  opacity 
values  to  pixels  based  on  features  defined  in  a  sep¬ 
arate  domain.  These  techniques  have  been  success¬ 
fully  applied  to  large  gigabyte-scale  data  sets  [8], 
[31].  Creating  useful  transfer  functions  for  visualizing 
three-dimensional  scalar  fields  is  understood  to  be  a 
difficult  problem.  In  these  cases,  the  primary  goal  is 
to  visualize  particular  structures  in  the  data  while 
minimizing  occlusion  from  surrounding  structures. 
The  most  common  technique  is  to  assign  color  and 
opacity  values  based  on  the  intensity  value  of  each 
pixel.  However,  complex  volumetric  data  sets  repre¬ 
sented  using  a  scalar  field,  such  as  magnetic  resonance 
imaging  (MRI)  and  computed  tomography  (CT),  often 
contain  multiple  tissue  types  at  similar  intensity  val¬ 
ues.  In  these  cases,  voxel  intensity  is  insufficient  for 
distinguishing  between  tissue  types.  Previous  work 
by  Kniss  et  al.  [19]  demonstrate  the  use  of  spatial 
characteristics  for  building  multi-dimensional  transfer 
functions.  The  transfer  function  domain  is  defined 
using  local  features,  such  as  pixel  intensity,  gradi¬ 
ent  magnitude,  and  curvature  [18].  Later  techniques 
demonstrate  the  use  of  transfer  functions  to  visualize 
volumetric  structures  based  on  large-scale  features, 
such  as  a  scale-space  representation  of  the  image  [7]. 
Techniques  have  also  been  developed  to  automate  the 
creation  of  transfer  functions  for  3D  images  based  on 


occlusion  and  the  spatial  distribution  of  data  [26],  [33]. 

However,  these  techniques  are  difficult  to  gener¬ 
alize  to  spectroscopic  images  since  each  pixel  repre¬ 
sents  a  spatially-resolved  absorbance  function.  Very 
few  techniques  currently  exist  for  visualizing  spec¬ 
troscopic  images.  Li  et  al.  [25]  propose  a  technique 
for  visualizing  astrophysical  data  imaged  at  vari¬ 
ous  wavelengths.  They  propose  the  use  of  transfer 
functions  for  defining  opacity  when  rendering  the 
image  stack  volumetrically.  However,  the  number  of 
bands  in  a  single  image  is  small  and  the  samples  are 
discontinuous,  which  limits  the  amount  of  chemical 
information  in  the  data  set.  More  recent  work  demon¬ 
strates  a  visualization  framework  for  near-infrared 
spectroscopic  images  of  historical  documents  sampled 
regularly  in  the  spectral  domain  [17].  This  technique 
allows  relighting  of  images,  interactive  selection  of  in¬ 
dividual  bands,  and  metrics  for  evaluating  similarity 
between  user-specified  spectra.  However,  similarity  is 
difficult  to  measure  for  vibrational  spectra  because 
of  the  coupled  structural  and  chemical  contributions 
to  the  spectrum.  Unsupervised  techniques,  such  as 
principal  component  analysis  (PCA)  and  vertex  com¬ 
ponent  analysis  (VCA)  [29],  have  also  been  proposed 
to  perform  spectral  unmixing.  However,  these  tech¬ 
niques  assume  homogeneous  samples  and  also  make 
specific  assumptions  about  the  data,  such  as  the  ex¬ 
istence  of  orthogonal  chemical  signatures,  which  are 
not  generally  applicable. 

Finally,  no  automated  techniques  have  been  pro¬ 
posed  that  compensate  for  spectral  features  intro¬ 
duced  by  elastic  scattering,  which  accounts  for  a 
large  amount  of  variance  in  mid-infrared  spectro¬ 
scopic  images  [22].  The  methods  that  we  propose  are 
independent  per-pixel  (per-spectrum)  operations  that 
allow  interactive  exploration  of  spectroscopic  images. 
We  use  two-dimensional  transfer  functions  based  on 
absorbance  and  wavelength,  combined  with  interac¬ 
tive  processing  techniques,  to  identify  and  visualize 
structural  and  chemical  information. 

4.1  Metrics 

We  first  identify  useful  measurements  for  quantifying 
spectral  features.  The  user  specifies  parameters  for 
these  measurements  in  the  TF  domain  using  a  simple 
user  interface.  A  color  map  is  then  applied  to  these 
metric  results  in  order  to  produce  a  transfer  function. 

4.1.1  Absorbance  at  a  Wavenumber  (Peak) 

The  absorbance  at  a  specified  wavenumber  is  the  most 
basic  spectral  measurement  and  is  given  by  A(x ,  y ,  v), 
where  (x,y)  is  the  spatial  position,  and  v  is  the 
wavenumber.  This  metric  is  extremely  fast  to  compute 
and  is  useful  when  the  chemistry  associated  with  a 
spectral  feature  is  well  understood.  It  is  also  highly 
sensitive  to  peak  shifts.  This  is  particularly  noticeable 
for  the  Amide  I  peak  at  1650cm-1,  which  is  narrow 
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Fig.  4.  Metrics  for  selecting  spectral  features.  Useful 
measurements  include  (a)  peak  height,  (b)  average 
peak  height  (integral),  (c)  centroid  (center  of  gravity). 
Colored  arrows  indicate  the  projection  of  the  metric 
value  on  the  widget  (black  line).  Shaded  regions  indi¬ 
cate  the  integral  intervals. 


and  composed  of  multiple  chemical  contributions  that 
make  it  prone  to  shift. 

Absorbance  is  particularly  useful  for  detecting  low- 
frequency  components,  like  scattering  and  broad 
peaks,  as  well  as  detecting  shifts  in  narrow  compo¬ 
nents.  However,  multiple  absorbance  values  must  be 
used  to  determine  the  direction  of  a  peak  shift.  It  is 
also  particularly  sensitive  to  noise. 

4. 1.2  Area  Under  the  Peak 

While  the  peak  absorbance  in  indicative  of  the  pres¬ 
ence  of  species,  the  total  area  under  the  peak  is 
actually  indicative  of  the  total  concentration.  Hence, 
another  useful  metric  is  the  definite  integral  of  the 
spectrum  on  a  user-specified  interval  [Po  ,  v\ }  that  pro¬ 
vides  the  area  under  the  peak: 

rv\ 

Mi(x,y,v o,Pi)=/  I(x,y,  v)dv  (5) 

J  D0 

This  metric  provides  a  robust  measure  of  absorbance 
within  a  spectral  region  and  is  relatively  insensitive 
to  noise  and  peak  shifts.  This  robustness  makes  it  an 
ideal  candidate  for  use  as  a  spectral  reference.  It  is 
insensitive,  however,  to  subtle  spectral  changes  that 
are  often  found  in  biological  tissue  samples.  For  color 
mapping  (Sec.  4.3),  the  average  peak  height  is  used 
instead  of  the  total  value  of  the  definite  integral. 


4.1.3  Centroid 


The  centroid  of  a  bounded  region  of  the  spectrum  is 
computed  using: 


Mc(x,y,v 0,1/1) 


fto  I(x,V,v)di> 


(6) 


The  resulting  value  is  the  wavenumber  for  the  center 
of  mass  in  the  specified  region.  This  metric  is  use¬ 
ful  for  measuring  shifts  in  single  peak  positions  as 
well  as  the  distribution  of  absorption  among  multiple 
neighboring  peaks.  This  metric  is  dependent  on  the 
distribution  of  absorption  and  therefore  does  not  re¬ 
quire  a  reference.  However,  it  is  incapable  of  detecting 
the  height  or  existence  of  peaks  that  do  not  shift  as  a 
function  of  spatial  position. 


4.2  Dynamic  Processing 

Our  proposed  TF  is  specified  in  a  space  defined 
by  absorbance  and  wavenumber.  The  data  set  is 
projected  into  this  domain  and  visualized  using  a 
joint  histogram.  The  user  can  then  select  features  in 
this  domain  that  correspond  to  the  structural  and 
chemical  components  of  the  data  set  to  be  visualized. 
However,  the  data  processing  required  to  separate 
coupled  structural  and  chemical  features  significantly 
affects  the  projection  of  the  data  into  the  TF  domain. 
We  implement  a  dynamic  projection,  which  allows 
the  user  to  explore  the  data  set  via  interactive  feed¬ 
back  in  both  the  spatial  and  spectral  domains,  which 
facilitates  meaningful  feature  selection.  Our  frame¬ 
work  allows  the  user  to  interactively  adjust  points  for 
baseline  correction  and  select  reference  features.  Data 
processing  is  performed  dynamically  on  the  GPU  us¬ 
ing  CUDA  [30]  and  provides  interactive  feedback  for 
multi-gigabyte  data  sets.  The  unprocessed  data  set  is 
stored  on  the  GPU  as  a  three-dimensional  texture  map 
represented  as  32-bit  floating  point  values.  The  user 
manipulates  the  data  in  two  ways:  (a)  the  insertion 
of  baseline  points  and  (b)  the  selection  of  a  reference 
metric. 

4.2. 1  Histogram 

The  joint  histogram  is  computed  by  assigning  a  sep¬ 
arate  block  of  threads  per  band  (wavenumber).  We 
specify  a  2D  block  size  of  ^/wxyjw  threads,  where 
w  is  the  maximum  warp  size  supported  by  the 
GPU.  Therefore,  each  block  consists  of  one  warp 
that  executes  data  in  a  single-instruction  multiple- 
data  (SIMD)  fashion.  Each  block  is  responsible  for 
computing  the  complete  one-dimensional  histogram 
for  a  single  band.  The  threads  within  each  block  are 
responsible  for  evaluating  a  spatially  coherent  square 
of  pixels,  initially  positioned  at  the  upper  left-hand 
corner  of  the  image.  Each  thread  then  iterates  across 
the  image  to  the  lower-right  corner  at  intervals  of  y7 w . 
Note  that  all  threads  within  the  block  are  part  of  the 
same  (SIMD)  warp,  therefore  they  will  be  spatially 
coherent  at  each  iteration  across  the  spatial  domain 
of  the  image. 

The  processed  spectrum  is  computed  dynamically 
as  the  user  changes  parameters  for  the  baseline  and 
reference.  Evaluation  of  the  processed  spectrum  re¬ 
quires  (a)  the  current  spatial  coordinates  and  band,  (b) 
the  bands  for  the  two  neighboring  baseline  points,  and 
(c)  the  reference  value.  The  band  is  provided  through 
the  block  identifier  (blockldx .  x),  while  the  current 
spatial  position  depends  on  the  2D  thread  identifier 
[x0,yo]  =  [threadldx .  x,  threadldx.y]  and  iter¬ 
ation  i,  where  [xo,yo\  is  the  initial  thread  position 
at  i  =  0.  The  baseline  bands  are  stored  in  constant 
memory  as  a  list.  Determining  the  two  neighboring 
baseline  bands  uq  and  v\  requires  a  binary  search, 
however  this  is  only  necessary  upon  block  initializa¬ 
tion  and  the  bands  are  stored  in  a  local  registers  for 
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Fig.  5.  Parallel  implementation  for  projecting  the  data  set  into  the  TF  space.  Each  band  is  evaluated  by  a  block 
of  y/wxy/w  threads,  (a)  Each  block  iterates  across  the  spatial  dimensions  of  the  image  from  i  =  0  to  i  =  n, 
where  n  =  and  Sx  and  Sy  are  the  spatial  extents.  The  threads  for  i  =  14  (*)  are  shown,  (b)  The  threads 
are  evaluated  using  SIMD  and  kept  spatially  coherent  to  take  advantage  of  2D  texture  caching,  (c)  Each  band  is 
evaluated  by  an  independent  thread  block  and  thread  performs  a  maximum  of  3  texture  fetches  per  iteration  (1 
data  point  and  2  baseline  points),  (d)  Resulting  histogram. 


all  iterations.  The  spectral  reference  is  pre-computed 
as  a  metric  (described  below)  and  stored  on  the  GPU 
as  a  2D  image. 

Computing  the  processed  spectrum  at  each  spatial 
location  within  a  band  requires  a  maximum  of  four 
memory  fetches:  (1)  the  raw  data  value  at  I(x,y,u), 
(2)  the  raw  values  at  the  baseline  points  I(x,y,v o) 
and  I{x,  y,  v\),  and  (3)  the  reference  value  r(x.  y).  This 
proposed  method  of  parallelization  provides  several 
advantages.  First  of  all,  the  use  of  independent  SIMD 
execution  for  each  band  allows  us  to  compute  the 
histogram  by  writing  only  to  shared  memory.  The 
histogram  is  then  copied  to  the  output  image  in  global 
memory  upon  termination  of  all  threads  in  the  block. 
By  using  only  one  warp  per  block,  we  also  force  all 
threads  to  perform  spatially  coherent  memory  fetches 
for  all  iterations.  This  allows  us  to  take  advantage  of 
spatially  coherent  2D  texture  caching  for  fetches  into 
both  the  raw  data  I  and  the  reference  image  r.  This 
is  an  efficient  use  of  GPU  resources,  provided  that 
the  number  of  bands  exceeds  the  number  of  available 
GPU  cores.  This  is  often  the  case,  however  idle  threads 
can  be  loaded  for  smaller  data  sets  by  subdividing 
the  spatial  domain  and  assigning  coherent  regions 
to  different  blocks.  This  will  require  performing  an 
atomic  merge  to  combine  the  resulting  histograms  in 
global  memory. 

Finally,  note  that  accurate  computation  of  the  his¬ 
togram  requires  the  use  of  atomic  adds  in  shared 
memory  in  order  to  prevent  conflicts  between  threads 
in  the  same  block.  However,  since  the  histogram  is 
used  for  display  purposes  only,  we  found  that  the 
use  of  atomic  operations  was  unnecessary  for  two 
reasons:  (1)  the  histogram  is  displayed  as  an  intensity 
map  and  the  display  device  only  has  an  8-bit  dynamic 
range,  (2)  intensity  values  where  conflicts  are  likely 
to  occur  are  dense,  making  them  less  noticeable  (Fig. 
6).  Therefore,  we  found  that  atomic  writes  to  shared 


(a)  (b)  (c) 

Fig.  6.  Comparison  of  histogram  rendering  (a)  with 
and  (b)  without  atomic  operations,  (c)  The  difference 
image  is  shown  (amplified  by  10X).  The  full  histogram 
is  shown  along  with  several  selected  spectra  (orange). 

memory  provided  little  improvement  in  quality  for 
large  data  sets  and  eliminating  them  gave  an  «4X 
increase  in  performance  (Sec.  5). 

4.2.2  Computing  Metrics 

Metrics  are  computed  independently  per  pixel  (per 
spectra),  therefore  spatial  partitioning  of  the  data  set 
does  not  introduce  the  need  for  atomic  operations. 
The  spectra  are  still  processed  dynamically,  requiring 
baseline  points  and  reference  bands.  The  data  set 
is  partitioned  spatially  into  y/wxy/w  blocks,  where 
each  thread  is  assigned  an  independent  spectrum. 
The  baseline  bands  and  reference  are  applied  as  de¬ 
scribed  previously  (Sec.  4.2.1).  If  the  metric  requires 
integration,  summation  occurs  along  the  spectral  axis. 
Note  that  a  binary  search  for  the  neighboring  baseline 
bands  is  only  necessary  when  a  thread  is  created. 
If  a  baseline  band  is  crossed  within  an  interval  of 
integration,  the  new  baseline  bands  are  determined 
by  acquiring  the  next  band  vn  in  the  list  and  shifting 
(To  Pi  P„).  The  metric  result  is  written  to  an 
image,  which  is  either  used  as  a  source  for  color¬ 
mapping  to  produce  the  final  display  image  (Sec.  4.3) 
or  as  a  reference  for  additional  metrics. 

4.3  Color  Mapping 

The  user  assigns  color  values  to  spectra  based  on 
measurements  using  a  simple  widget.  Color  mapping 
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(b)  (c) 


Fig.  7.  Color  mapping  of  spectral  features  in  a 
baseline-corrected  image.  Labels  indicate  peaks  rep¬ 
resenting  prominent  chemical  features,  (a)  The  spec¬ 
tral  histogram  shows  three  intervals  of  constant  color 
(RGB)  applied  at  the  widget  position  at  Amide  I  (vertical 
white  line).  Note  how  the  colors  mix  as  individual 
spectra  cross  and  overlap  with  changes  in  density 
and  chemistry.  Spatial  images  show  the  result  of  a  (b) 
Gaussian  and  (c)  linear  ramp  applied  to  the  widget. 

is  incorporated  into  the  selection  of  metric  parameters. 
The  user  places  the  widget  at  the  desired  metric  posi¬ 
tion  in  the  histogram.  For  all  three  metrics,  the  widget 
position  is  indicated  by  a  line  segment.  For  metrics 
that  are  computed  across  a  specified  wavenumber 
interval,  a  shaded  region  indicating  these  bounds  is 
also  displayed  (Fig.  4). 

When  the  result  of  a  spectral  measurement  inter¬ 
sects  this  line  segment,  a  color  value  is  assigned  based 
on  the  point  of  intersection  using  a  ramp  shader. 
Our  framework  allows  the  user  to  assign  colors  using 
a  constant,  linear,  and  Gaussian  ramp  (Fig.  7).  In 
order  to  make  selection  intuitive  for  the  integration 
metric  M,  ,  color  is  assigned  based  on  the  average  peak 
height: 


This  allows  the  widget  to  be  placed  in  the  vicinity  of 
the  selected  spectra. 

5  Results 

Rendering  time  is  dominated  by  the  dynamic  com¬ 
putation  of  the  histogram  as  the  user  changes  pro¬ 
cessing  and  visualization  parameters.  However,  the 
frame  rate  is  interactive  for  our  largest  sample  image 
(700x700x491,  «1GB),  requiring  <  32ms  for  complete 
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(a)  (b) 


(c)  (d) 


Fig.  8.  Basic  demonstration  of  the  proposed  pro¬ 
cessing  and  visualization  methods  on  a  USAF  test 
target  constructed  from  SU-8  photoresist  on  a  BaF2 
substrate,  (a)  Magnitude  of  the  absorbance  vector  and 
the  full  histogram  (logarithmic  intensity  scaled),  (b) 
Scattered  spectra  selected  and  displayed,  (c)  Baseline 
corrected  spectra  showing  polymer  thickness  using  a 
linear  map.  Note  the  region  of  dual-layered  polymer 
(arrow),  (d)  Normalized  spectra  using  the  metric  from 
(b)  as  a  reference.  Features  associated  with  SU-8  are 
selected. 


evaluation  of  the  histogram  (<  147ms  with  atomic 
writes  to  shared  memory).  Testing  was  performed  on 
a  GeForce  GTX  580  with  1.5GB. 

In  this  section,  we  use  the  proposed  techniques  to 
develop  transfer  functions  for  for  both  mid-infrared 
and  Raman  spectroscopic  imaging  data.  We  first 
demonstrate  the  visualization  of  structural  and  chem¬ 
ical  components  from  images  of  well-characterized 
polymer  targets  that  are  often  used  to  assess  image 
quality.  We  then  show  how  these  techniques  can  be 
used  to  extract  similar  information  from  mid-infrared 
images  of  tissue  biopsies,  including  the  visualization 
of  chemical  components  useful  for  breast  cancer  diag¬ 
nosis.  Finally,  we  show  how  the  proposed  techniques 
can  be  used  on  Raman  images  to  visualize  the  distri¬ 
bution  of  proteins  and  sugars  in  soybean  images. 
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Fig.  9.  Using  transfer  functions  to  visualize  features  seen  in  standard  histology.  A  surgical  resection  of  a  breast 
biopsy  was  imaged  using  IR  spectroscopy  and  then  stained  and  imaged  with  a  standard  optical  microscope,  (a) 
A  section  of  the  hematoxylin  and  eosin  (H&E)  stained  tissue  sample  is  shown.  The  proposed  techniques  can  be 
used  to  selectively  visualize  various  structural  and  chemical  features  using  only  the  IR  hyperspectral  image,  (b) 
Scattering  effects  are  highlighted,  indicating  regions  containing  highly  scattering  structures  (edges  and  collagen 
fibers),  (c)  Removal  of  scattering  features  allows  the  selection  of  tissue  thickness  based  on  the  absorbance  of 
Amide  I  (1650cm-1).  (d)  Using  Amide  I  as  a  reference,  chemical  compounds  corresponding  to  epithelial  cells 
can  then  be  visualized,  (e-f)  The  corresponding  joint  histograms  are  shown,  with  selected  pixels  highlighted. 


5.1  Multi-Component  Air  Force  Target 

We  demonstrate  the  basic  principles  behind  the  pro¬ 
posed  transfer  functions  using  an  image  of  a  known 
sample.  We  constructed  a  1951  USAF  resolution  test 
chart  using  an  SU-8  photoresistive  polymer  on  a  bar¬ 
ium  flouride  (BaF2)  substrate.  The  target  is  acquired 
using  a  commercial  mid-infrared  spectroscopic  imag¬ 
ing  system  (Perkin-Elmer  Spotlight  300)  and  visual¬ 
ized  using  a  colormap  based  on  vector  magnitude  and 
compared  to  images  constructed  using  the  proposed 
transfer  functions  (Fig.  8).  We  are  able  to  cleanly 
extract  pixels  associated  with  edges  from  the  raw  data 
due  to  enhanced  scattering  in  their  spectra.  After  base¬ 
line  correction,  we  can  then  estimate  polymer  thick¬ 
ness  if  the  absorption  coefficients  are  known  or  use 
a  normalization  peak  to  examine  relative  absorbance 
for  the  general  case.  Note  that  the  absorbance  of 
pixels  at  the  edges  decreases  as  expected,  while  the 
a  double-layered  region  exhibits  approximately  twice 
the  absorbance.  We  then  select  a  reference  peak  at 
1120cm-1  (a  component  of  SU-8)  for  normalization 
and  select  another  SU-8  feature  in  order  to  extract 
the  pixels  that  compose  the  target  image.  All  spectral 
histograms  are  shown  using  log-scale  intensity. 

5.2  Infrared  Images  of  Tissue 

Biological  tissue  samples  pose  a  unique  problem  in 
spectroscopic  imaging,  since  each  pixel  contains  a 
complex  combination  of  chemical  species.  Spectral 
features  that  correspond  to  different  tissue  types  are 
subtle  and  difficult  to  identify,  making  interactive 
exploration  a  valuable  tool.  We  demonstrate  the  use  of 


the  proposed  methods  for  the  identification  of  tissue 
types  in  mid-infrared  images. 

We  first  demonstrate  the  structural  and  chemical 
characteristics  that  can  be  visualized  using  this  tech¬ 
nique.  We  obtain  a  tissue  sample  and  place  it  on  an 
IR-reflective  glass  substrate  for  imaging.  The  same 
sample  is  then  stained  using  a  commonly-used  clin¬ 
ical  dye.  Hematoxylin  and  Eosin  (H&E),  in  order  to 
identify  structural  and  chemical  features  (Figure  9a). 

The  proposed  method  is  used  to  mine  similar  fea¬ 
tures  from  the  hyperspectral  image  data.  Since  the 
unprocessed  spectra  are  dominated  by  scattering  in 
the  high-wavenumber  regions,  we  are  able  to  visu¬ 
alize  pixels  containing  edges  and  boundaries,  where 
scattering  is  more  prominent  (Fig.  9b  and  e).  We  then 
select  baseline  points  and  use  a  linear  gradient  to 
visualize  the  tissue  thickness  and  density  based  on 
Amide  I  absorption  (Fig.  9c  and  f),  in  which  over  97% 
of  cancers  of  the  breast  arise.  The  facile  delineation 
of  epithelial  cells  is  a  desirable,  but  unmet  need,  as 
reported  in  recent  studies  using  automated  comput¬ 
erized  analysis  of  pathologic  images  [4]. 

Finally,  we  demonstrate  the  ability  to  identify  ep¬ 
ithelial  cells  in  the  sample,  which  exhibit  unique 
chemical  characteristics.  By  using  Amide  I  as  a  ref¬ 
erence,  features  in  the  spectra  are  now  dominated  by 
chemical  differences  between  cell  types  in  the  tissue. 
This  allows  the  user  to  visualize  the  distribution  of 
tissue  types,  such  as  epithelial  cells  (Fig.  9d  and  g). 

We  further  demonstrate  the  possible  use  of  the  pre¬ 
sented  transfer  functions  for  tissue  pathology.  Several 
breast  cancer  biopsies  were  commercially  obtained 
and  stained  using  H&E.  Adjacent  tissue  samples  were 
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(d)  (f)  (g) 


Fig.  1 0.  Comparison  of  tissue  visualization  techniques,  (a)  Standard  histology  slide  stained  with  H&E.  Epithelial 
cells  are  purple,  connective  and  necrotic  tissue  is  stained  pink,  (b)  Spectral  magnitude  indicated  using  a  rainbow 
color  map.  (c)  The  first  three  principle  components  displayed  embedded  in  the  RGB  color  channels,  (d)  Chemical 
composition  and  density  using  spectral  TFs,  and  (e)  the  corresponding  spectral  histogram.  Widget  position 
indicates  the  source  of  (d),  and  the  shaded  region  is  used  as  a  reference,  (f)  Referenced  histogram  and 
corresponding  (g)  chemical  visualization  of  the  tissue  sample.  Features  correspond  to  collagen  (red),  epithelium 
(green),  and  necrosis  (blue).  Arrows  indicate  corresponding  features  in  the  histology.  Scale  bar  represents 
100/im.  The  IR  data  set  is  700x700  pixels  with  490  spectral  samples  per  pixel. 


imaged  using  mid-infrared  spectroscopy.  We  then  de¬ 
sign  transfer  functions  to  identify  tissue  characteristics 
that  are  useful  for  cancer  diagnosis.  This  includes 
structural  features,  such  as  tissue  density,  as  well  as 
labeling  of  various  tissue  types.  These  tissue  types 
include  epithelial  cells,  connective  tissue  (stroma),  and 
necrotic  (dead  and  dying)  cells.  The  spatial  distribu¬ 
tion  of  these  features  are  commonly  used  in  pathology 
to  identify  the  presence  of  a  tumor.  Finally,  we  com¬ 
pare  our  proposed  methods  to  those  commonly  used 
for  spectroscopic  image  analysis,  including  unsuper¬ 
vised  techniques  such  as  PC  A  (Fig.  10). 

The  ability  to  identify  tissue  types  by  interactively 
exploring  the  hyperspectral  data  set  provides  several 
advantages  over  existing  techniques.  Unlike  chem¬ 
ically  staining  the  tissue  to  extract  knowledge  of 
constituent  histology,  mid-infrared  imaging  does  not 
perturb  the  tissue  sample  in  any  way.  In  addition. 


the  effects  of  elastic  scattering  on  tissue  samples 
are  not  well  understood,  and  previous  methods  for 
classification  of  tissue  types  require  extensive  pre¬ 
processing  and  the  application  of  machine  learning 
[14],  [27].  Though  an  extensive  comparison  between 
the  more  complex  machine  learning  methods  and  the 
method  reported  here  must  be  undertaken  to  quan¬ 
tify  advantages,  the  results  seem  to  be  comparable 
and  significantly  faster  for  our  method  by  facilitat¬ 
ing  interactive  analysis.  In  contrast,  an  experienced 
user  was  able  to  separate,  visualize,  and  evaluate  the 
utility  of  the  extracted  tissue  characteristics  useful  for 
tumor  diagnosis  in  a  few  minutes,  given  interactive 
feedback. 

5.3  Raman  Spectroscopy 

Some  advantages  of  Raman  spectroscopy  include  the 
ability  to  choose  the  excitation  frequency,  easy  cou- 
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(C) 


Fig.  11.  Raman  spectroscopy  of  a  soybean,  (a)  A 
standard  rainbow  colormap  of  spectral  magnitude,  (b) 
Visualization  showing  the  distribution  of  protein  (red), 
carbohydrates  and  unsaturated  fats  (green)  and  lipids 
(blue),  (c)  The  Raman  histogram  is  shown  with  promi¬ 
nent  peaks  labeled. 

pling  to  fiber-optic  probes,  and  little  to  no  interference 
by  the  existenc  of  water  in  a  sample.  These  features 
make  it  useful  for  taking  in  vivo  measurements  [36], 
[35],  [34].  However,  the  Raman  spectral  signal  is  ex¬ 
tremely  weak  and  requires  long  acquisition  times  that 
effectively  limit  the  signal-to-noise  ratio  (SNR)  and  the 
size  of  image  that  can  be  obtained  in  a  reasonable 
amount  of  time.  In  this  section,  we  demonstrate  that 
transfer  functions  can  be  created  to  effectively  visual¬ 
ize  the  distribution  of  chemical  compounds  in  Raman 
images.  We  collected  spectra  from  a  cross-section  of 
a  soybean  and  visualized  the  distribution  of  spectral 
bands  correlated  with  sugars,  proteins,  and  oils  (Fig. 
11). 

6  Conclusions  and  Future  Work 

The  ability  to  acquire  spatially  resolved  information 
using  hyperspectral  imaging  is  strongly  emerging  as  a 
promising  avenue  across  a  variety  of  areas,  especially 
biomedical  analyses.  Vibrational  spectroscopy  has  the 
potential  to  provide  quantitative  histology  for  disease 
diagnosis,  without  the  use  of  chemical  stains  in  an 
objective  and  automated  manner.  However,  computa¬ 
tional  methods  must  be  developed  to  visualize  the 
data  quickly  and  reliably  The  size  and  complexity 
of  the  data  contained  in  hyperspectral  images  makes 
this  a  difficult  problem,  requiring  the  separation  of 
physical  and  chemical  characteristics  from  underlying 
spectra.  In  this  paper,  we  demonstrate  an  interactive 
method  for  building  transfer  functions  for  visualizing 
hyperspectral  images.  Our  method  allows  users  to 
dynamically  assimilate  large  collections  of  spectra 
using  algorithms  designed  to  separate  structural  and 
chemical  features  in  real  time.  These  features  are 
then  selected  using  transfer  functions  which  allow 
the  visualization  of  these  characteristics  in  a  spatial 
image  of  the  sample.  To  our  knowledge,  the  reported 
methods  are  the  first  to  allow  interactive  processing 


and  visualization  of  hyperspectral  images  at  this  level 
of  spectral  and  structural  detail,  and  we  have  demon¬ 
strated  their  usefulness  in  biological  samples. 

Future  directions  include  applying  these  techniques 
to  three-dimensional  samples,  which  can  be  acquired 
using  Raman  spectroscopy  in  combination  with  a  con- 
focal  microscope.  However,  acquiring  these  images  is 
currently  extremely  time  consuming  due  to  the  weak 
signal  obtained  using  Raman  imaging  methods. 

Our  proposed  technique  also  has  several  advan¬ 
tages  over  unsupervised  methods,  such  as  PCA  and 
VCA.  In  particular,  the  metrics  that  we  use  for  vi¬ 
sualization  have  finite  support,  requiring  only  a  nar¬ 
row  band  of  information  within  the  spectrum.  Once 
useful  metrics  are  identified,  the  number  of  collected 
bands  can  then  be  reduced,  allowing  faster  imaging, 
for  example  using  narrow-band  filters  for  IR  [20]. 
Finally,  since  the  separation  of  structural  and  chemical 
characteristics  from  an  IR  image  is  so  difficult,  many 
algorithms  for  the  classification  of  hyperspectral  im¬ 
ages  rely  on  the  use  of  user-defined  metrics  [14].  Our 
method  may  provide  an  efficient  method  for  selecting 
features  for  use  in  more  complex  classifiers  for  IR- 
based  clinical  histology. 
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Dear  Editor, 

We  are  submitting  our  manuscript  titled  "Designing  Transfer  Functions  for  Interactive  Data  Mining  of 
Hyperspectral  Images"  for  consideration  for  publication  in  IEEE  Transactions  on  Visualization  and 
Computer  Graphics. 

This  work  describes  an  interactive  method  for  exploring  hyperspectral  images,  which  are  emerging 
in  several  fields  including  astronomy,  chemistry,  and  biomedicine.  We  place  a  particular  emphasis 
on  immediate  applications  for  cancer  diagnosis  in  biomedical  images. 

Please  note  that  a  version  of  this  work  was  accepted  as  a  poster  presentation  at  IEEE  VisWeek 
2011.  We  have  attached  the  corresponding  2-page  paper.  All  of  the  text  and  images  have  been 
updated,  and  the  paper  provides  significantly  more  background  information.  We  have  included 
necessary  implementation  details  as  well  as  a  thorough  description  of  our  technique  and 
applications  on  several  additional  data  sets. 

Thank  you  for  your  consideration, 

David  Mayerich,  Michael  Walsh,  Matt  Schulmerich,  Rohit  Bhargava 
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Designing  Transfer  Functions  for  Exploring  Hyperspectral  Images 

David  Mayerich  Michael  Walsh  Rohit  Bhargava 

Beckman  Institute  for  Advanced  Science  and  Technology,  University  of  Illinois  at  Urbana-Champaign 


Abstract — Spectroscopy  is  used  in  several  fields  to  acquire  chemical  information  about  a  substance  for  identification  or  analysis. 
Recent  advances  allow  the  acquisition  of  hyperspectral  imagery  for  biomedical  data.  This  has  promising  implications  in  the  fields  of 
biotechnology  and  medicine,  where  quantitative  analysis  of  tissue  can  be  performed  by  directly  measuring  spatially-resolved  chemical 
information.  However,  this  requires  a  significant  amount  of  human  intervention  to  identify  spectral  features  for  use  in  classification.  We 
propose  a  method  for  designing  transfer  functions  for  hyperspectral  images.  This  method  allows  researchers  to  interactively  adjust 
parameters  used  to  manipulate  the  input  spectra  in  order  to  find  metrics  that  can  be  used  to  classify  features  in  the  images. 

Index  Terms — biomedical  imaging,  medical  imaging,  spectroscopy,  transfer  functions 

-  ♦  - 


1  Introduction 

Hyperspectral  data  is  composed  of  a  series  of  samples,  generally  taken 
across  the  electromagnetic  spectrum.  Different  compounds  exhibit  a 
characteristic  signature  in  the  spectral  vector,  which  can  be  used  to 
determine  the  chemical  properties  of  a  sample.  This  technique  is  of¬ 
ten  used  in  the  fields  of  biomedicine  and  forensics  to  determine  the 
identity  of  chemical  compounds. 

Recent  advances  in  detector  technology  allow  hyperspectral  signals 
to  be  resolved  spatially,  producing  multidimensional  images  with  an 
additional  spectral  component.  These  techniques  have  shown  promise 
in  biomedicine  where  sections  of  tissue  can  be  classified  by  cell  type 
based  on  the  spatially  resolved  spectra  from  a  mid-infrared  image  [2]. 
However,  the  interpretation  of  a  spectral  signature  is  a  complex  task 
that  currently  requires  human  intervention  to  identify  spectral  features 
for  use  in  classification.  In  addition,  a  significant  amount  of  pre¬ 
processing  must  be  applied  before  identifying  spectral  components 
that  correspond  to  structural  and  chemical  information.  This  is  be¬ 
cause  chemical  and  structural  features  of  the  tissue  are  not  represented 
as  independent  spectral  components.  Therefore,  determining  chemi¬ 
cal  features  often  requires  the  removal  of  structural  properties,  such  as 
cell  density  and  scattering  characteristics. 

The  goal  of  this  work  is  to  create  a  tool  to  aid  spectroscopists  in 
selecting  features  for  tissue  classification.  Current  automated  meth¬ 
ods,  such  as  principle  component  analysis  (PCA)  and  vertex  compo¬ 
nent  analysis  (VCA)  [5]  operate  on  variance  and  are  prone  to  cap¬ 
turing  noisy  features,  such  as  spectra  distorted  through  scattering.  In 
addition,  the  components  produced  using  PCA  and  VCA  have  broad 
support  across  the  entire  spectrum.  Since  IR  imaging  can  be  time- 
consuming,  finding  localized  metrics  is  important  for  clinical  applica¬ 
tions,  since  they  allow  classification  with  fewer  spectral  samples.  Ac¬ 
curate  classification  methods  combined  with  fast  and  non-destructive 
chemical  imaging  can  play  an  important  role  in  clinical  research,  par¬ 
ticularly  for  the  diagnosis  of  cancer  biopsies. 

We  present  a  method  that  allows  interactive  exploration  of  hyper¬ 
spectral  data  sets  by  building  transfer  functions  similar  to  those  used 
for  volume  rendering.  We  demonstrate  preliminary  results  using  bi¬ 
ological  data  acquired  from  mid-infrared  spectroscopic  imaging.  We 
first  define  a  transfer  function  showing  the  distribution  of  spectra  in  the 
data  set.  We  then  allow  the  user  to  specify  color  values  in  this  domain. 


•  David  Mayerich  is  with  the  Beckman  Institute,  University  of  Illinois  at 
U rbana- Champaign,  mayerich  @  illinois.  edu. 
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•  Rohil  Bhargava  is  with  the  Department  of  Bioengineering  and  the 
Beckman  Institute,  University  of  Illinois  at  Urbana-Champaign. 


Fig.  1 .  Transformations  applied  to  hyperspectral  data  in  order  to  extract 
chemical  information,  (left)  Baseline  estimation  (green)  removes  arti¬ 
facts  due  to  light  scattering.  Normalization  to  known  peaks,  such  as 
Amide-1  (red)  compensates  for  tissue  density  and  thickness,  (a)  Distri¬ 
bution  of  spectra  that  have  been  baseline  corrected  and  normalized  to 
Amide-1  and  (b)  the  2D  image  of  the  Amide-1  peak  before  normalization 
(after  baseline  correction).  A  threshold  value  is  specified  at  the  normal¬ 
ization  wavelength  in  order  to  prevent  the  introduction  of  artifacts  from 
division  by  small  numbers. 


This  technique  is  similar  to  those  published  previously  using  spatial 
characteristics  [4].  However,  we  also  allow  the  user  to  dynamically 
adjust  preprocessing  parameters,  thereby  changing  the  distribution  of 
features  in  the  spectral  domain. 

1.1  Mid-Infrared  Spectroscopy 

Mid-infrared  (IR)  spectroscopic  imaging  is  performed  by  measuring 
the  amount  of  light  absorbed  by  a  tissue  sample  at  wavelengths  in  the 
mid-IR  region,  generally  between  2 jim  and  12/im.  The  characteris¬ 
tic  signature  of  a  chemical  compound  is  caused  by  the  absorption  of 
photons  by  molecular  bonds  found  in  organic  molecules.  The  level 
of  absorption  at  any  point  in  the  spectrum  can  be  correlated  with  the 
existence  of  certain  atomic  bonds.  These  correlations  can  be  used  to 
identify  the  compound  being  analyzed.  Other  features,  such  as  the 
scattering  properties  of  the  tissue  sample  and  the  cell  density  at  a  given 
spatial  location,  also  affect  the  spectrum  and  may  provide  additional 
features  useful  for  visualization. 

1.2  Image  Analysis 

The  absorbance  spectrum  measured  at  a  particular  spatial  location  in 
a  tissue  sample  is  affected  by  three  major  factors:  (a)  scattering  of 
light  through  the  tissue  sample,  (b)  the  tissue  thickness  and  density, 
and  (c)  the  chemical  composition  of  the  material.  Scattering  effects 
are  caused  by  the  tissue  sample  acting  as  an  imperfect  lens,  distorting 
the  light  as  it  passes  through  the  specimen.  Since  scattering  is  based 
on  the  wavelength  of  transmitted  light,  this  distortion  differs  along  the 
spectrum.  In  general,  this  results  in  a  shift  in  the  baseline  and  peak 
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(d)  (e)  (f) 


Fig.  2.  Regions  highlighted  in  a  log-scale  joint  histogram  of  the  spec¬ 
tral  domain  with  various  preprocessing  steps  applied,  (a,  d)  Scattering 
effects  allow  edges  to  be  highlighted  by  selecting  spectra  in  the  raw 
data  in  a  region  known  to  contain  no  chemical  information,  (b,  e)  Af¬ 
ter  a  baseline  is  applied,  cell  density  is  highlighted  by  setting  a  gradient 
along  the  O-H  bonding  region,  (c,  f)  Normalization  of  the  spectra  to  the 
Amide  I  peak  allows  selection  of  more  complex  chemical  features,  such 
as  DNA-rich  epithelial  cells. 


positions.  While  mathematical  models  have  been  proposed  to  com¬ 
pensate  for  scattering  effects  [1],  they  rely  on  an  understanding  of  the 
expected  output  spectrum.  For  unknown  tissue  samples,  a  piecewise 
linear  baseline  correction  is  generally  applied  (Figure  1). 

Tissue  density  and  thickness  cause  scaling  of  the  absorbance  ac¬ 
cording  to  Beer’s  Law: 

Ax=b'£euxci  (1) 

i 

where  A  is  the  absorbance  measured  at  a  spatial  location,  b  is  the  path 
length  (tissue  thickness),  q  is  the  concentration  of  compound  i  and 
£i  ^  is  the  absorbance  of  compound  i  at  wavenumber  X  per  unit  of 
concentration  [3].  Since  Amide  I  (found  at  «  1650cm-1)  is  a  protein 
present  in  most  tissue  samples,  spectra  are  often  normalized  to  the 
value  at  this  wavenumber  for  initial  analysis.  However,  peak  height 
ratios  are  often  useful  metrics,  therefore  normalization  to  other  peaks 
for  exploration  is  important. 

2  Transfer  Function  Design 

We  display  the  projection  of  the  data  set  into  the  transfer  function 
domain  using  a  two-dimensional  joint  histogram  of  spectral  distribu¬ 
tion  as  a  function  of  spectral  component  X  and  signal  amplitude.  The 
histogram  for  each  slice  is  computed  on  the  graphics  processor  us¬ 
ing  CUDA.  Each  value  from  the  raw  data  set  is  transformed  based 
on  user-specified  parameters  for  baseline  correction  and  normaliza¬ 
tion.  The  projection  of  the  data  set  into  the  transfer  function  domain  is 
computed  interactively,  allowing  the  user  to  dynamically  change  pre¬ 
processing  parameters  and  view  the  corresponding  changes  in  the  dis¬ 
tribution  of  spectra.  This  allows  the  user  to  segregate  chemical  and 
structural  features  of  the  tissue  into  localized  regions  in  the  spectral 
domain.  The  user  can  then  select  these  features,  which  are  rendered  in 
a  two-dimensional  orthographic  projection  of  the  data  set  in  the  spa¬ 
tial  domain  (Figure  2  a-c).  The  distribution  of  selected  spectra  are  also 
shown  as  an  overlay  in  the  spectral-domain  image  (Figure  2  d-f). 

A  typical  search  for  classification  metrics  involves  labeling  known 
compounds  in  the  spatial  domain.  This  can  be  done  using  histology 
slides  of  neighboring  sections  as  a  reference.  The  user  then  tunes  the 
baseline  and  normalization  criteria  to  find  parameters  that  localize  the 
selected  spectra  at  some  point  in  the  spectral  domain.  The  user  tests 
this  region  as  a  metric  by  highlighting  it  using  a  series  of  widgets  [4] 
and  validating  the  results  in  the  2D  spatial-domain  projection. 

3  Results 

We  demonstrate  preliminary  results  of  this  technique  by  using  dy¬ 
namic  hyperspectral  transfer  functions  to  segment  physical  and  chem¬ 


Fig.  3.  (a)  A  tissue  sample  imaged  using  IR  and  stained  with  hema¬ 
toxylin  and  eosin  (H&E)  for  histology.  The  H&E  section  is  compared  to 
the  same  region  visualized  with  a  transfer  function  specified  using  the 
proposed  method.  Tissue  features  such  as  cell  density  (red),  edge  ef¬ 
fects  (green)  and  epithelium  (blue)  are  selected  independently  and  com¬ 
bined  into  a  2D  visualization  (b)  and  in  the  spectral  domain  (c).  While  the 
selection  of  cell  types,  such  as  epithelium,  require  baseline  correction 
and  normalization,  these  preprocessing  steps  eliminate  other  features 
from  the  spectra  as  shown  by  the  overlapping  regions  in  the  transfer 
function  (c). 


ical  characteristics  of  data  from  a  breast  biopsy.  We  show  edge  pixels 
prone  to  scattering  artifacts  by  selecting  them  from  the  unprocessed 
spectra.  By  applying  baseline  correction,  we  determine  tissue  density 
based  on  IR  absorption  at  the  Amide  I  band.  By  normalizing  to  Amide 
I,  we  then  locate  a  spectral  region  unique  to  epithelial  cells  that  line 
lumen  in  the  tissue  sample.  The  segmented  cells  and  corresponding 
spectral  regions  are  shown  in  Figure  2  and  the  final  classified  image 
and  spectral  overlays  are  shown  in  Figure  3.  As  seen  in  the  full  spec¬ 
tral  overlay  (Figure  2b),  independent  selection  of  these  components 
would  be  impossible  without  manipulating  the  preprocessing  parame¬ 
ters. 

Because  of  the  size  and  complexity  of  these  data  sets,  the  ability  to 
interactively  search  multiple  parameter  spaces  provides  a  much  more 
efficient  method  of  metric  selection.  For  future  work,  we  will  test  met¬ 
rics  selected  using  this  technique  against  those  computed  using  other 
methods,  such  as  PC  A,  YCA,  and  purely  manual  selection.  In  addition 
to  providing  metrics  localized  to  small  regions  of  the  spectrum,  we  ex¬ 
pect  that  the  number  of  required  metrics  can  be  reduced,  improving 
classification  by  limiting  over-fitting. 
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Abstract 

Image  quality  from  an  infrared  microscope  has  traditionally  been  limited  by  considerations  of  throughput 
and  signal  to  noise  ratio  (SNR).  A  first  principles  understanding  of  the  achievable  quality  as  a  function  of 
instrument  parameters,  however,  is  needed  for  improved  instrument  design.  Here,  we  first  present  a  model 
for  light  propagation  through  an  infrared  (IR)  spectroscopic  imaging  system  based  on  scalar  wave  theory. 
The  model  analytically  describes  the  propagation  of  light  along  the  entire  beam  path  between  the  source 
and  the  detector.  The  effect  of  various  optical  elements  and  the  sample  in  the  microscope  is  understood  in 
terms  of  the  accessible  spatial  frequencies  using  a  Fourier  optics  approach  and  simulations  are  conducted  to 
gain  insights  into  spectroscopic  image  formation.  The  optimal  pixel  size  at  the  sample  plane  is  calculated 
and  shown  to  be  much  smaller  than  that  in  current  mid-infrared  microscopy  systems.  A  commercial  imaging 
system  is  modified  and  experimental  data  are  presented  to  demonstrate  the  validity  of  the  developed  model. 
Building  on  this  validated  theoretical  foundation,  an  optimal  sampling  configuration  is  set  up.  Acquired  data 
were  of  high  spatial  quality  but,  as  expected,  of  poorer  SNR.  Signal  processing  approaches  are  implemented 
to  improve  the  spectral  SNR  and  resulting  data  demonstrate  the  ability  to  perform  high- definition  IR  imaging 
in  the  laboratory  using  minimally- modified  commercial  instruments. 
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1  Introduction 


Instrumentation  for  Fourier  transform  infrared  (FT-IR)  microspectroscopic  imaging  typically  consists  of  an 
interferometer  for  multiplexed  spectral  encoding,  microscope  optics  for  condensing  light  and  image  formation 
as  well  as  a  focal  plane  array  (FPA)  detector  for  multichannel  data  recording. 1  The  instrumentation  has 
benefited  from  nearly  60  years  of  development  with  its  genesis  in  point-by-point  mapping,2,3  and  FT-IR 
microscopy  developed  in  the  1980s.4  The  first  instruments  contained  a  single-element  detector  that  collected 
ah  light  transmitted  by  the  microscope  while  apertures  were  used  to  define  the  spatial  resolution.  The  use 
of  far-held  apertures  to  localize  the  region  illuminated  at  the  focal  plane  of  the  microscope  implied  that  the 
smallest  spot  size  attained  was  primarily  determined  by  the  wavelength  of  light.  Most  studies,  however, 
involved  larger  spot  sizes  to  achieve  higher  throughput  and,  consequently,  higher  signal  to  noise  ratio  (SNR) 
of  acquired  data.  Two  dogmatic  ideas  then  emerged  to  dominate  IR  imaging  technology.  The  first  was  that 
significant  information  could  not  be  derived  from  areas  smaller  in  dimension  than  the  wavelength  of  light. 
While  this  was  indeed  true  due  to  a  lack  of  throughput  for  the  point  microscopy  case,  the  elimination  of 
apertures  in  IR  imaging  meant  that  there  would  still  be  significant  throughput  at  smaller  pixels  sizes.  The 
second  misconception  was  that  there  was  no  benefit  to  increasing  pixel  density  beyond  the  optical  diffraction 
limit  imposed  by  the  wavelength  as  it  would  not  result  in  more  spatial  details  being  accessible.  Evidence  to 
the  contrary  was  available  with  the  use  of  sub-wavelength  apertures  along  with  the  higher  throughput  of  a 
synchrotron  source,  making  the  approach  feasible  and  providing  excellent  quality  data. 5 

In  other  studies,  there  was  also  evidence  of  improved  image  quality  by  sub-pixel  stepping6  in  map¬ 
ping  experiments  that  did  not  gain  much  favor  on  two  grounds.  The  first  objection  was  a  practical 
consideration.  Long  acquisition  time  became  even  longer.  The  second  objection  arose  due  to  lack  of  an 
appropriate  theory,  leading  to  a  mixing  of  the  concepts  of  the  resolution  of  optical  microscopy  with  the 
quality  of  images  arising  from  IR  absorption.  Often,  the  two  are  used  interchangeably,  which  we  emphasize 
in  this  manuscript,  is  not  the  case.  Though  it  is  correct  that  a  resolution  higher  than  determined  by  the 
optical  configuration  and  wavelength  cannot  be  attained,  there  is  also  no  denying  the  improvement  in  image 
quality  by  recording  far-held  data  from  areas  of  dimensions  smaller  than  wavelength.  For  example, 7-9 
aperture  sizes  as  small  as  3 fim  x  3 fim  provided  for  excellent  data  using  a  synchrotron  source  at  wavelengths 
much  longer  than  3  iim.  Spectral  quality  was  retained  in  this  study  while  maps  were  sharper.  More 
recent  studies10  have  shown  improved  image  quality  by  combining  sub-pixel  images.  However,  the  data 
themselves  have  not  been  collected  at  the  optimal  pixel  sizes.  Coupling  of  an  FPA  to  synchrotron-based 
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microscope  systems6,11,12  has  recently  provided  stunning  improvements  in  image  quality. 13  Spectra  of 
reasonable  quality  could  be  acquired  in  minutes  by  using  multiple  beams  and  the  multichannel  advantages 
of  FPAs  allowed  for  large  area  coverage.  In  these  advances,  the  brightness  of  a  synchrotron  has  been  a 
key  factor  and  it  is  unclear  if  the  same  benefits  may  be  achievable  in  conventional  instruments  using  a 
globar  source.  While  there  have  been  considerable  advances  in  instrumentation,  few  theoretical  analyses 
and  validation  of  resolution  and  imaging  properties  have  been  undertaken. 14,15  These  analyses,  however, 
did  not  provide  for  a  rigorous  model  of  the  IR  microscope  and  did  not  examine  the  case  of  pixel  sampling 
beyond  the  conventional  wavelength-limited  designs.  The  lack  of  an  appropriate  theory  to  optimize  image 
quality,  hence,  remains  as  a  gap  in  our  understanding.  This  gap  is  also  a  barrier  to  design  of  table-top 
instruments  of  optimal  quality  as  it  is  difficult,  at  present,  to  make  design  trade-offs  in  a  quantitative  manner. 

In  this  manuscript,  we  first  undertake  a  theoretical  analysis  of  the  image  formation  in  IR  microspec- 
troscopic  imaging.  In  particular,  we  focus  on  determining  the  best  configuration  for  optimal  image  quality. 
We  emphasize  that  attainable  resolution  and  image  quality,  though  related  by  wavelength  and  parameters  of 
the  optical  setup,  may  not  be  determined  by  the  same  rules.  Once  the  optimal  image  quality  is  determined, 
only  then  can  an  examination  of  the  resolution  attainable  be  undertaken.  Hence,  a  detailed  discussion  of 
resolution  is  not  undertaken  here  and  we  focus  instead  on  criteria  for  obtaining  images  of  the  highest  possible 
quality.  Using  the  developed  theory,  we  first  determine  the  pixel  size  for  highest  image  quality  that  is 
permitted  by  the  optical  components.  While  we  focus  here  only  on  the  capabilities  of  the  microscope,  there 
are  other  theoretical  approaches  dealing  comprehensively  with  the  effects  of  optical  properties  of  the  sample, 
size  and  shape  of  domains,16  geometry  of  acquisition  and  microscopy  optics  on  both  spectral  distortions17,18 
and  image  formation. 19  The  approach  here,  hence,  is  a  more  general  model  that  neglects  sample  detail  in 
the  interests  of  understanding  microscope  performance.  Second,  we  implement  the  recommendations  from 
theory  by  modifying  a  commercial  system.  Finally,  we  present  an  integration  with  previously  described 
noise-reduction  strategies  to  provide  improved  data  quality. 

2  A  Model  for  FT-IR  Spectroscopic  Imaging 

A  theoretical  analysis  is  presented  here  for  the  system  shown  in  Fig.  1.  Models  based  on  ray  (geometric) 
optics  have  been  especially  popular  for  IR  microscopy,  but  are  inadequate  for  analysis  of  wavelength 
associated  effects.  Full-scale  electromagnetic  models,  at  the  other  extreme,  capture  all  of  the  physics  of 
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image  formation  including  polarization  effects,  but  also  involve  increased  complexity.  It  is  prudent  to  choose 
the  simplest  model  that  incorporates  all  the  phenomena  of  interest.  Hence,  we  forgo  the  full  electromagnetic 
treatment  in  favor  of  using  the  simpler  scalar  wave  theory  framework. 

Our  analysis  is  based  on  explicitly  calculating  the  electric  field  at  every  plane  in  the  microscope  sys¬ 
tem.  To  simplify  the  analysis,  we  consider  a  monochromatic  component  of  the  field  with  a  wavenumber  v 
and  complex  amplitude  U .  Dependence  on  v  is  implied  throughout  unless  otherwise  stated.  The  field  may  be 
expressed  as  a  superposition  of  plane  waves  described  by  functions  of  the  form  U  (f)  exp  {i2ir  [f  •  r  +  fz(f)z]} 
where  f  is  a  two-dimensional  vector,  for  instance  in  Cartesian  coordinates  f  =  (fx,fy)-  The  function  fz( f) 
is  defined  so  that  |f|2  +  /2( f)  =  z/2,  that  is,  the  plane  waves  satisfy  the  Helmholtz  equation.  For  each 
plane  wave  constituting  the  field,  its  relative  contribution,  determined  by  its  amplitude,  and  its  relative 
position,  determined  by  its  phase  are  both  necessary  for  a  complete  description  of  the  field.  The  field  and 
all  linear  transformations  of  the  field  can  be  represented  as  a  linear  combination  of  these  plane  waves,  which 
travel  at  various  angles.  The  term  “angular  spectrum”  is  often  used20  to  describe  such  a  decomposition. 
We  emphasize  that  the  use  of  the  word  “spectrum”  in  this  context  should  not  be  confused  with  the 
absorption  spectrum  that  is  commonly  encountered  in  spectroscopy  and  that  “spectrum”  without  the 
modifier  “angular”  is  meant  to  be  the  usual  absorption  spectrum.  Similarly,  the  term  “spatial  frequency” 
refers  to  the  vector  f  used  above  and  should  not  be  confused  with  spectral  “frequency” .  For  the  convenience 
of  the  reader,  we  have  summarized  the  notation  used  in  appendix  A. 2. 

2.1  Optical  System 

The  system  shown  in  Fig.  1  can  be  analyzed  using  a  modular  approach.  To  that  end,  we  describe  each 
component  in  the  optical  system  with  an  operator.  This  approach  provides  a  convenient  means  of  modifying 
the  analysis  to  accommodate  changes  in  the  instrument.  Since  the  system  is  linear,  it  is  convenient  to  work 
in  a  notation  designed  for  linear  algebra,  namely  the  Dirac  notation.21  Vectors  describing  the  field  in  a  plane 
at  constant  axial  coordinate  2  are  written  | Uz)  with  Hermitian  adjoint  (Uz\.  The  two-dimensional  Fourier 
component  of  the  field  at  spatial  frequency  f  is  written  U( f)  =  (f\Uz)  and  the  value  of  the  field  at  a  point 
r  is  given  by  U( r)  =  (r\Uz).  The  sample  is  assumed  to  be  thin,  as  described  by  a  vector  | S).  The  detector 
plane  is  taken  to  be  at  z&.  Thus  the  goal  of  the  analysis  is  to  produce  an  equation  of  the  form 

\UzD)  =  A  (S') ,  (1) 
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Fig.  1.  Model  for  a  FT-IR  imaging  system.  The  system  consists  of  an  interferometer  coupled 
to  a  microscope  system.  Cassegrain  C\  focus  light  on  to  the  sample  S',  C2  and  Cd  transmit 
light  from  the  sample  to  the  detector  D.  The  fields  in  different  planes  are  also  indicated. 

where  A  describes  the  optical  system. 

We  first  address  the  operator  for  propagation  through  free  space  between  two  parallel  planes  sepa¬ 
rated  by  a  distance  d,  denoted  IQ.  We  assume  unidirectional  propagation  along  the  axis  normal  to  the 
planes  so  that  propagation  is  simply  the  accumulation  of  phase.  The  matrix  elements  of  IQ  are  given  by 


<fi|  Kd  |f2)  =  <S(f2  -  h)ei2nMfl)d, 


(2) 


where  S  is  the  Dirac  delta  function.  Alternatively,  the  operator  may  be  constructed 


(3) 


It  is  useful  to  note  that  IQ  IQ/  =  K d+d'  • 


An  interferometer  may  be  thought  of  as  a  beam  splitter  followed  by  propagation  of  two  different  ax¬ 
ial  distances.  Given  arms  of  length  d\  and  measured  from  the  beam  splitter  to  each  of  the  mirrors,  the 


interferometer  is  represented  by  an  operator 


Idi ,d2  =  ^  (K2di  +  K 2d2)  ■ 


(4) 


Next  we  consider  the  effect  of  a  Cassegrain.  Much  like  a  Fresnel  lens,  the  Cassegrain  can  be  considered  to 
impart  onto  the  field  a  quadratic  phase  factor.  We  represent  the  Cassegrain  C  with  an  operator  Gc, 


J  d2r  |r)  Qc(r)exP 


.iris  n 
-i—r 

Lf 


(rl 


(5) 


where  Lf  is  the  focal  length  of  the  Cassegrain.  The  quadratic  phase  factor  converts  each  plane  wave  incident 
on  the  Cassegrain  into  a  spherical  wave  converging  on  its  focal  plane.  The  phase  relation  between  plane 
waves  incident  on  the  Cassegrain  along  with  the  imparted  quadratic  phase  determine  the  final  image  position 
and  size.  Qc  is  the  aperture  function  (described  in  detail  in  the  appendix).  It  is  often  the  case  that  the 
Cassegrain  appears  between  two  propagation  steps  in  the  form  K^G^K^.  As  shown  in  the  appendix,  to 
a  good  approximation,  when  d^1  +  d^1  =  TJ1,  the  Cassegrain  focuses  light  from  the  first  plane  into  the 
second  with  a  magnification  factor  Me  =  cfe/di*  Thus,  it  is  convenient  to  define 


Hc  =  Kd2GcKdl  »  J  d2f  |f)  Qc  (f)  (~MC f  | . 


(6) 


Clearly  He  is  not  uniquely  defined  unless  d\  and  d 2  are  specified  and  some  care  should  be  taken  in 
observing  the  context  in  which  it  is  placed.  When  the  image  plane  is  out-of- focus  by  a  distance  do,  for 
example,  the  propagation  operator  K^0  may  be  used  to  obtain  the  out-of-focus  image  using  the  operator 
Kd0+d2GcK^  =  Kdollc. 


Given  a  field  produced  at  the  source  |C/o) ,  at  2  =  zo  then  propagated  to  the  interferometer  at  z  =  zj, 
acted  on  by  the  interferometer,  then  propagated  to  the  first  Cassegrain,  Ci,  at  2  =  zqx  and  then  from  the 
Cassegrain  to  the  focal  plane  of  the  Cassegrain,  z  =  Zf±,  the  resultant  field  is 

|  UZfl)  =  ^Zf1-ZClQc1KZc1-zI^-d1,d2^zi-ZQ  |f^o)  (7) 

The  sample  is  modeled  as  a  thin  screen.  That  is,  the  field  on  the  far  side  of  the  sample  is  assumed  to  be  the 
product  of  the  sample  structure  function  and  the  field  incident  on  the  sample  in  the  coordinate  domain.  We 
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formally  construct  an  operator  that  effects  this  multiplication.  Given  a  field  |  [/),  the  operator  U  is  defined 


u  =  I d2r  |r)  <r| U)  (r|  =  J d 2/d2/'  |f)  (f  -  t'\U)  <f'| .  (8) 

The  sample,  |S),  then  interacts  with  the  incident  field,  | UZfi ) ,  to  produce  a  transmitted  field  at  a  plane 
z  =  Zf1+  just  past  the  sample 


UZ}1+)  =  TJZ/1  \S)  =  J  d2/ d2/'  |f>  (f  -  f'l^)  (f'|S> . 


(9) 


Substituting  from  Eq.  7  we  find  that 

\UZfi  +  )  =  J d2/d2/'  |f)(f-f/|KZ/i_ZCiGc1K,Ci_Z/Idl)d2K2,_,0|C/o)(f'|5).  (10) 

Note  that  Id1:d2  =  ^2d1^-o,d2-d1-  Assuming  no  misalignment,  i.e.  when  the  system  is  in  focus,  the  propa¬ 
gation  operators  can  be  absorbed  into  the  magnification  factor  of  H c±  •  Thus 

| UZfi+)  =  Jd2fd2f  |f)(f-f'|HClIM2_dl|[70)(f'|S)  (11) 

where  =  ~KZf  zc  Gc1’KZCl-z0+2d1-  Evaluating  Eq.  11  explicitly,  we  find  that  the  field  just  past  the 
sample  is  given  by 

| UZfi+)  =  ^Jd2fd2f\f)Qc1(f'/Mc1){l+exp\i2nfz(fl)(2d2-2d1)}}(f/\Uo)(f  +  f,/Mc1\S)  (12) 

where  it  may  be  noted  that  (f  \S)  is  the  Fourier  component  of  the  sample  at  spatial  frequency  f.  In  the  case 
that  the  source  consists  of  a  point  source  far  from  the  first  Cassegrain,  (f  \Uo)  ~  Ao,  where  inclination  factors 
of  the  form  1  /  fz  have  been  neglected.  Then  the  field  is  given  by 

| UZJ1  +  )  =  ^  J  d2fd2fl\f)Qc1(f'/MCl){l+exp[i27rfz(i,)(2d2-2d1)}}(i  +  i,/MCl\S).  (13) 

Light  propagation  from  the  Cassegrain  C\  to  the  sample  and  to  the  Cassegrain  C2  are  illustrated  in  Fig.  2. 
The  propagation  of  the  field  after  the  sample,  \UZf  +},  to  the  field  at  the  detector,  | UZD ) ,  takes  place  through 
two  more  Cassegrains.  The  Cassegrain  C2,  located  effectively  at  z  =  zc2 ,  focuses  the  sample  plane  to  a  plane 
z  =  zf2  such  that  ^  \zc — b  where  Lf2  is  the  focal  length  of  Cassegrain  C2.  Similarly,  the 


7 


Fig.  2.  Cassegrain  Ci,  which  focuses  light  on  to  the  sample,  and  C2,  the  collection 
Cassegrain,  are  shown. 


Cassegrain  Cd  located  at  z  —  zcD ,  focuses  light  from  the  plane  z  =  ZfD  onto  the  detector  plane  at  z  =  Zd 

such  that  ^  x_z^ - h  Zd  _}Zc  =  z~~-  Misalignment  of  the  focal  planes  such  that  ZfD  7^  Zf2  may  be  taken 

into  account  with  the  propagator  -Zf2  •  Thus 


\UZD)  =  UcDKZfD_ZfHc2\UZfi  +  ). 

=  I d2/d2/'  HCDKZfD_ZfHc2  |f)  (f  -  f'|  HClI0,d2_dl  \U0)  (f'|S> .  (14) 

This  expression  is  now  straight-forward  to  evaluate.  We  assume  that  ZfD  =  Zf2  (the  focal  planes  are  aligned). 
We  find  explicitly  that 


I UZD)  =  \jd2fd2f'\i)QcDmcA-MCDf)QCl{-{'/MCl) 

x  {1  +  exp[i2n  f z(f')  (2d?  -  2d,)}}  (f'\U0)  (f'/MCl  +  MC2MCJ\S) .  (15) 

Define  M  =  McDMc2  and  compute 

(Wzo)  =  \QcD(f)Qc2(-McDf)  I d2fQc1(-i,/MCl) 

x  {1  +  exp[i27r/,(f,)(2d2  -  2d!)]}  (i'\U0)  (?/ MCl  +  Mf\S) .  (16) 
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The  above  expression  is  true  for  any  spatially  heterogeneous  sample  and  has  no  restriction  on  its  constituent 
spatial  frequencies.  Note  that  M  is  typically  larger  than  1  (i.e.  the  detector  elements  are  bigger  than  the 
corresponding  areas  on  the  sample). 


In  order  to  illustrate  the  significance  of  Eq.  16  consider  a  case  where  the  sample  is  a  point  object. 
In  this  case  (f '  /Mc1  +  Mf  \S)  =  1,  giving 

<f| UZD)  =  l-QcD{f)QcA-McDf)  J dPf'QcA-t'/Mc,) 

x  {1  +  exp [i2n fz (f ') (2d2  -  2dx)]>  <f'|E7„>  (17) 

=  \QcD(t)Qc2(-McDt)B0.  (18) 

Here,  Ho,  the  integral  from  Eq.  17  is  independent  of  f.  The  corresponding  detector  intensity  is 

(UZD\f)  (i\UZD)  =  ^[QcD({)Qc2(-MCDf)]*[QCD(i)QC2(-MCDf)}  (19) 

where  *  represents  convolution. 

In  a  case  where  the  source  is  completely  incoherent,  i.e.  (  (UZo  |r/)  (r\UZo)  )e  =  Io(r)5(r  —  r'),  the 
detector  intensity  is 


(  (UZD  |r)  (r| UZD)  )e  =  J  d2rl0(r')\h(r-r')\2  (20) 

where  h( r;  r')  is  the  transfer  function  from  the  coherent  case  and  /o(r)  is  the  source  intensity.  The  subscript 
e  above  denotes  an  ensemble  average  over  constituent  random  processes.  Using  the  equations  presented 
in  this  section,  we  proceed  to  compute  the  spatial  sampling  rate  and  pixel  size  required  to  record  all  the 
information  available  from  the  FT-IR  spectroscopic  imaging  system. 


2.2  Pixel  Size 

An  analysis  of  the  equations  for  the  FT-IR  imaging  system  provides  insight  into  optical  design  and  image  for¬ 
mation.  The  explicit  form  of  Qct  in  terms  of  NA  is  derived  in  appendix  A.l  (Eqs.  29-31).  It  may  be  observed 
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from  Eq.  16  that  the  spatial  frequencies  on  the  detector  are  always  less  than  min(u'N Ac2 out Ac DOUt) 
because  the  support  of  the  data  is  the  intersection  of  the  support  of  Qc2(~McD f)  and  QcD{ f)-  When  the 
pre-optics  to  the  detector  are  well  designed,  then  z/NA c2out/M  <  z/NA cDout  providing  the  limiting  field  on 
the  detector  to  be  a  spatial  frequency  of  z/NA c2out/M.  Since  the  detector  records  intensity  and  not  fields,20 
the  intensity  /( r)  =  E*(r)E(r)  has  a  spatial  frequency  bound  of  2  x  z/NA c2out/M.  In  order  to  faithfully 
record  the  entire  intensity  image  without  any  loss  of  information,  the  spatial  sampling  rate  according  to 
Nyquist  criterion20  has  to  be  at  least  twice  this  limiting  frequency,  or  4z/NA c2out/M.  The  pixel  size  is 
inversely  related  to  the  sampling  rate,  i.e.  Lpixei  =  M/(49NAc20Ut)-  The  wavenumber,  z/,  in  mid-infrared 
spectroscopic  imaging  typically  varies  between  600  cm-1  and  4000  cm-1  (0.4  Therefore,  using  Eq. 

16,  the  equivalent  pixel  size  on  the  sample  has  to  be  1/(4  x  0.4  NA^  out)  /tm.  Thus,  for  a  NA c2out  =  0-5, 
the  effective  pixel  size  on  the  sample  is  1.25  /rm  and  for  NA c2out  =  0.65,  the  effective  pixel  size  is  0.96  /im. 
It  must  be  emphasized  that  the  value  of  pixel  size  calculated  above  is  not  the  same  as  resolution.  It  is 
in  fact  smaller  than  the  resolution  as  defined  by  the  Rayleigh  criterion.20  However,  it  reflects  the  least 
sampling  rate  required  to  utilize  all  the  information  passed  by  the  system  in  FT-IR  spectroscopic  imaging. 

While  we  have  focused  on  image  quality,  a  number  of  other  insights  are  also  apparent.  In  some  de¬ 
tectors,  for  example,  it  may  be  advantageous  to  have  a  concave  mirror  instead  of  a  Cassegrain  and  that 
can  be  easily  be  incorporated  into  the  model  by  making  the  term  corresponding  to  the  central  obscuration 
i.e.  inner  numerical  aperture  zero  NA^-^0.  The  absence  of  a  central  obscuration  increases  light  throughput 
significantly.  We  suggest  that  the  Cassegrain  used  in  commercial  instruments  be  replaced  with  this  mirror. 
This  modification  will  result  both  in  better  image  quality  and  lower  instrument  costs.  Although  the  above 
derivation  uses  a  transmission  mode  measurement,  the  same  analysis  can  be  performed  for  transflection 
mode.  In  the  case  of  transflection-absorption  measurement  using  a  beam-splitter  and  full  Cassegrain  illumi¬ 
nation,  Cassegrains  C\  and  C2  correspond  to  the  same  Cassegrain  and  light  travels  through  the  sample  twice, 
once  before  reflection  and  once  after.  However,  the  sampling  rates  and  pixel  sizes  are  the  same  in  both  cases. 


3  Experiments  and  Simulations 

Instrumentation.  A  Varian  7000  Spectrometer  coupled  to  UMA-400  microscope  was  used  to  perform 
two  sets  of  experiments  with  parameters  listed  in  Table  1.  System  1  uses  accessories  standard  to  the 
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Varian  microscope.  System  2  uses  an  Edmund  Optics  NA=0.65,  74 x  magnification  Cassegrain  for  C2. 
The  two  systems  make  comparisons  easy  since  parameters  other  than  Cassegrain  C2  are  the  same.  The 
only  difference  is  in  the  effective  pixel  size  on  the  sample  and  NA c2out-  The  interferometer  is  operated 
in  the  step-scan  mode  at  a  stepping  rate  of  200  Hz.  Data  were  acquired  at  every  other  zero-crossing  (a 
undersampling  ratio  of  2)  of  a  He-Ne  laser  for  a  free-scanning  spectral  range  of  7900  —  0  cm-1.  A  Fourier 
transform  of  the  recorded  data  was  carried  out  using  the  Norton-Beer  medium  apodization  function.  Data 
were  truncated  and  stored  as  absorbance  after  a  ratio  against  an  appropriate  background  data  set. 


Parameter 

System  1 

System  2 

NA  clOUt 

0.5 

0.5 

NA  c2ont 

0.5 

0.65 

NA  Coout 

0.5 

0.5 

MCi 

15x 

15x 

MC2 

15x 

74  x 

mCd 

3x 

3x 

Pixel  Size  (effective) 

5.5/rm 

1.115  /rm 

Spectral  Resolution 

8cm  _1 

8cm  _1 

#  Scans  per  pixel 

8 

128 

Scans  (Background) 

128 

128 

Detector  size 

128  x  128 

128  x  128 

Undersampling  ratio 

2 

2 

Table  1.  Experimental  Parameters  used  for  data  acquisition  and  modeling. 


Data  extraction  and  image  processing  were  performed  using  a  hyper-spectral  imaging  software  package, 
ENvironment  for  Visualizing  Images  (ENVI).  In  ENVI,  the  two  dimensional  Fourier  transform  routines  are 
in  built.  Modules  for  further  data  processing  were  written  in-house  with  the  use  of  IDL  7.1  and  Matlab 
7.  The  two  imaging  systems,  with  parameters  shown  in  table  1  (parameters  1-6)  are  simulated.  The 
simulations  were  performed  on  a  quad  core  computer  with  a  Nvidia  GeForce  GTX  580  GPU  and  12GB 
RAM.  Modules  and  functions  for  the  simulations  were  written  in  house  using  Matlab  7.  The  simulations 
compute  the  detector  field  and  intensity  for  incoherent  illumination  using  the  system  parameters  (table  1, 
parameters  1-6)  and  object  as  inputs.  The  run  time  for  simulations  is  about  23s  for  each  wavenumber. 

Prostate  and  breast  tissue  microarrays  (TMAs)  were  used  as  a  platform  for  high  throughput  sam¬ 
pling. 22  TMAs  consist  of  a  large  number  of  small  (=  1mm  diameter)  tissue  sections  arranged  in  a  grid 
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pattern  on  an  substrate.  This  technique  facilitates  and  streamlines  acquisition  of  data  from  several  patients 
and  provides  diversity  to  the  sample  population  under  observation.  Breast  and  prostate  tissue  microarrays 
were  obtained  from  the  University  of  Illinois  Chicago  (kindly  provided  by  Dr  Andre  Kajdacsy-Balla)  and 
from  Biomax  (No.  BR1003;  US  Biomax,  Rockville,  MD).  One  5  fim  thick  tissue  section  from  each  array 
was  placed  on  a  BaF2  substrate  for  IR  imaging  and  serial  5  fim  thick  sections  were  placed  on  standard  glass 
slides  for  hematoxylin  and  eosin  ( H&E )  staining.  IR  imaging  data  were  acquired  from  a  normal  breast 
tissue  core  that  contained  a  region  comprising  of  terminal  ductal  lobular  units  as  well  as  from  a  normal 
prostate  tissue  core  from  the  region  of  a  small  blood  vessel. 


4  Results  and  Discussions 

Our  first  step  in  analyzing  the  performance  of  the  imaging  systems  was  to  validate  simulations  with 
recorded  data  using  a  gold  standard  for  comparison.  As  common  in  other  forms  of  microscopy,  a  standard 
USAF  1951  target,  consisting  of  chrome  on  glass,  was  chosen.  The  chrome  provides  a  large  attenuation 
(absorbance)  with  negligible  thickness.  It  also  provides  high  contrast  for  visible-microscope  imaging.  Hence, 
it  is  an  ideal  sample  for  simulation  and  validation.  We  first  recorded  visible-light  images  for  a  standard 
USAF  microscope  target  (Edmund  optics).  The  recorded  image  was  binarized  to  remove  any  edge  blurring 
or  grayscale  values  in  the  simulation.  An  absorbance  of  unity  was  assigned  to  the  chrome  region  and  the 
resulting  image  was  used  for  simulation.  Regions  without  the  chrome  were  assigned  a  transmittance  of 
unity  (zero  absorbance).  Simulations  were  carried  out  for  both  the  instrument  configurations  described 
previously.  The  light  source  is  assumed  to  be  spatially  incoherent  as  in  Eq.  20.  Spectroscopic  imaging 
data  were  also  recorded  using  the  target.  We  analyzed  group  number  6  and  7,  elements  1  through  6,  as 
these  represent  the  highest  image  quality  typically  encountered  in  IR  microscopy.  Images  of  the  standard 
sample  employed  for  simulation,  absorbance  images  from  simulated  data  and  absorbance  images  from 
the  recorded  data  corresponding  to  the  two  instrument  configurations  are  shown  in  Fig.  3.  In  order  to 
facilitate  comparison  of  simulations,  we  have  used  the  same  pixel  size  for  both  optical  configurations, 
although  different  NA  lenses  are  often  associated  with  different  magnifications.  The  chosen  pixel  size  for 
simulations  of  0.36  fim  will  not  be  a  limiting  factor  in  the  resulting  image  quality.  The  USAF  target  is 
designed  to  enable  facile  calculation  of  the  frequency  response  of  the  imaging  systems  as  well  as  to  illustrate 
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resolution  capabilities  via  images.  In  simulations,  we  emphasize  again  that  the  sample  is  assumed  to  be 
infinitesimally  thin  to  avoid  effects  associated  with  thickness.  Comparing  the  two  simulated  absorbance 
images  (Fig.  3(b)  and  (f))  with  the  same  pixel  size  shows  the  effects  of  increasing  NA.  As  expected,  an 
increase  in  the  NA  provides  for  higher  quality  and  higher  resolution  images.  Nevertheless,  the  enhancement 
is  not  especially  striking.  A  more  obvious  degradation  in  image  quality  is  evident  when  the  pixel  size  is 
increased  from  1.1  fim  (Fig.  3(f))  to  5.5  /im  (Fig.  3(g))  while  all  other  parameters  are  held  constant. 
This  dramatic  change  in  image  quality  emphasizes  the  need  for  carefully  choosing  the  pixel  size.  The 
combined  effects  of  choosing  a  higher  NA  lens  and  the  appropriate  sampling  (small  pixel  size)  can  result 
in  significantly  improved  image  quality  (cf.  Fig.  3(b)  and  Fig.  3(g)).  To  experimentally  validate  these 
predicted  improvements  in  image  quality,  we  set  up  two  optical  configurations  described  in  table  1  on  the 
same  commercial  microscope  and  obtained  data  on  the  same  USAF  1951  target  used  in  simulations.  Our 
goal  was  to  minimally  modify  (by  changing  one  lens)  the  existing  commercial  system.  Hence,  we  did  not 
alter  the  condenser.  Matching  the  throughput  via  the  image  formation  lens  and  condenser  is  likely  to 
result  in  better  SNR  but  will  not  materially  change  the  appearance  of  the  images  as  the  image  quality 
depends  only  on  the  angular  spectral  bandwidth  of  the  image  formation  lens.  Experimental  results  for  (c) 
NA=0.65,  pixel  size  =  1.115  //m  and  (h)  NA=0.50,  pixel  size  =  5.5  /rra  are  presented  for  comparison.  A 
dramatic  improvement  in  image  quality  is  observed  experimentally  for  the  smaller  pixel  size.  Expanding  a 
smaller  region  for  comparison  in  Figs.  3(d), (e)  and  3(i),(j),  the  three  smaller  bars  are  distinguishable  only 
in  the  higher  NA,  smaller  pixel  size  setup.  This  qualitatively  demonstrates  the  gain  in  image  quality  and 
resolution  using  a  higher  NA  system  along  with  the  appropriate  pixel  spacing.  In  all  cases,  as  expected, 
we  also  note  that  the  recorded  data  are  of  somewhat  lower  quality  than  the  simulations.  This  is  due  to 
edge  effects  in  the  sample,  finite  detector  sizes,  imperfect  optics  and  experimental  noise  in  the  recorded  data. 

While  the  images  from  simulation  and  recorded  data  agree  well,  they  do  not  provide  a  quantitative 
understanding  of  the  performance  improvement  and  limits  of  the  two  configurations.  Therefore,  we 
quantitatively  analyzed  the  performance  of  the  systems  as  a  function  of  the  spatial  frequency.  The  frequency 
response  of  the  optical  system  is  often  characterized  by  the  modulation  or  contrast  transfer  function 
of  the  system.  While  optical  systems  typically  employ  intensity  values  to  examine  the  imaging  system 
performance,  here  we  use  the  spatial  frequencies  in  Fig.  3(a)  and  the  corresponding  absorbance  values  from 
simulations  and  recorded  data  to  plot  the  absorbance  contrast.  The  absorption  contrast  ratio  (ACR)  is 
defined  as  the  difference  in  maximum  and  minimum  absorbance  observed  at  a  specific  spatial  frequency 
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compared  to  that  observed  in  the  ideal  case.  The  ideal  difference  between  high  and  low  absorbance  is 
obtained  by  comparing  absorbance  values  from  regions  on  the  sample  with  no  absorbance  (for  example, 
the  region  used  to  collect  a  background  spectrum)  and  a  relatively  homogeneous  region  of  high  absorbance 
(for  example,  on  the  large  square  in  a  USAF  1951  target).  Figure  4(a)  shows  a  plot  of  the  absorbance 
contrast  ratio  as  a  function  of  spatial  frequency  at  3950  cm-1.  As  expected,  ACR  is  highest  at  low  spatial 
frequencies  and  decreases  with  increasing  spatial  frequency  until  zero  (no  measurable  contrast).  ACR 
for  the  lower  NA  system  reaches  zero  at  a  lower  spatial  frequency  than  that  for  the  high  NA  one  as  is 
evident  from  the  polynomial  fit  (solid  lines)  to  simulation  data.  It  may  be  seen  in  the  experimental  data 
that  the  high  NA  system  resolves  all  bars  including  the  smallest  bars  (Group  #7)  available  on  the  USAF 
1951  target.  The  frequency  at  which  ACR  reaches  zero  can  be  predicted  via  simulations  to  be  about  0.4 
/mi-1  (=400  line  pairs/mm).  Contrast  decreases  with  reduction  in  pixel  size  since  light  intensity  per  pixel 
reduces.  As  a  result,  the  increased  imaging  capability  is  partially  offset  by  a  decreased  spectral  SNR  in 
the  resulting  absorbance  images.  While  predictions  from  simulations  agree  well  with  measured  ACR  for 
the  high  NA  system,  there  is  some  disagreement  for  the  lower  NA  setup  at  high  spatial  frequencies.  This 
difference  is  observed  at  approximately  0.08  fim -1  (=80  line  pairs/mm).  Experimentally,  the  lower  NA 
system  cannot  resolve  bars  beyond  a  spatial  frequency  of  0.115  iim~x  (=115  line  pairs/mm).  This  is  lower 
than  the  theoretically  predicted  0.15  fim~x  (=150  line  pairs/mm).  This  can  be  attributed  to  the  finite  pixel 
size  preventing  an  accurate  recording  of  the  data.  This  comparison  between  ACR  of  the  high  and  low  NA 
systems,  enabled  by  the  developed  theoretical  understanding,  provides  convincing  proof  that  the  designs  of 
present-day  instrumentation  are  not  permitting  optics-limited  performance  and  can  be  rapidly  and  easily 
reconfigured  to  provide  significantly  higher  imaging  performance. 


Though  the  previous  analysis  provides  the  smallest  observable  features,  the  relationship  of  these  fea¬ 
tures  to  the  optimal  pixel  size  needs  to  be  quantified.  We  examine  the  pixel  size  at  the  sample  plane 
required  for  optimally  sampling  the  signal  as  a  function  of  NA  and  wavelength  (Fig.  4(b)).  While  the  optimal 
pixel  size  can  also  be  calculated  analytically,  this  plot  provides  the  design  parameters  for  instruments 
should  a  practitioner  develop  an  instrument  for  high  performance  at  any  specific  wavelength  while  using  a 
specific  lens.  As  expected,  the  recommended  pixel  size  decreases  with  decreasing  design  wavelength.  The 
dependence  on  NA  is  rather  striking  as  well.  It  is  noteworthy  that  most  commercially  available  systems 
provide  reasonably  high  NA,  but  the  pixel  size  is  not  commensurate  with  the  optical  components.  Thus, 
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Fig.  4.  (a)  A  plot  of  absorbance  contrast  ratio  (ACR)  measured  from  both  simulation  and 
recorded  experimental  data  as  a  function  of  spatial  frequency  at  v  =  3950  cm-1  for  the 
two  systems  in  table  1.  (b)  A  plot  of  minimum  pixel  size  required  for  correctly  sampling 
allowable  spatial  frequencies  as  a  function  of  v  and  NA  is  shown.  The  pixel  size  corresponds 
to  the  effective  pixel  size  at  the  sample. 


an  improvement  in  incorporating  a  higher  NA  lens  is  not  likely  to  be  as  dramatic  as  using  the  optimal 
pixel  size  in  current  commercial  systems.  We  note  that  this  plot  is  based  on  the  maximum  attainable 
spatial  frequencies  for  thin  objects  and  larger  sample  thicknesses  are  likely  to  achieve  poorer  contrast  for 
the  limiting  spatial  frequency. 

The  pixel  sizes  recommended  will  provide  the  highest  image  quality  provided  by  the  microscope,  un¬ 
less  limited  by  the  sample.  The  optimal  pixel  size  and  the  resolution  of  an  imaging  system  are  also  related. 
According  to  the  Rayleigh  criterion,  for  example,  the  attainable  resolution  is  0.61  /(pNA).  We  propose  to 
sample  at  0.25/(z /NA)  in  order  to  convert  the  analog  signal  into  a  digital  readout  without  information  loss. 
We  note  that  the  pixel  size  proposed  is  significantly  smaller  than  the  resolution.  This  result  is  contrary  to 
prevailing  wisdom  in  that  it  is  widely  believed  that  no  further  improvement  is  possible  once  a  pixel  size  equal 
to  the  resolution  has  been  attained.  Indeed,  mixing  the  concepts  of  resolution  and  pixel  sizes  for  correct 
signal  sampling  has  led  to  significant  confusion  in  IR  microscopy.  We  emphasize  that  our  suggestions  for 
improved  data  quality  are  not  a  contravention  of  the  resolution  criterion.  Consider  the  case  of  pixel  spacing 
resulting  in  optimal  resolution.  To  resolve  two  objects,  we  need  more  than  two  measurements  when  the  data 
are  digitized.  For  example,  consider  two  point  objects  that  are  centered  on  a  pixel  each.  To  entirely  resolve 
these  objects,  appropriate  sampling  implies  that  there  be  a  pixel  on  either  side  of  the  object  to  separate 
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Fig.  5.  Absorbance  images  from  spectroscopic  imaging  data  obtained  from  the  two  config¬ 
urations  (NA=0.65,  NA=0.50)  are  shown  on  top.  The  corresponding  angular  spectra  are 
shown  at  the  bottom.  Images  for  the  system  with  (a)  NA=0.65  and  (b)  NA=0.50  obtained 
by  plotting  the  absorbance  at  v  =  2962  cm-1.  Similarly,  the  absorbance  at  v  =  1650  cm-1 
is  shown  for  systems  with  (c)  NA=0.65  data  and  (d)  NA=0.50. 


them  from  each  other  and  any  other  neighboring  objects.  Therefore,  at  least  5  pixels  are  needed  to  resolve 
two  point  objects.  A  more  formal  discussion  of  sampling  for  optical  microscopy6  and  its  extension  to  IR 
microscopy  is  provided  elsewhere  using  this  signal  processing  approach.  Our  approach  is  distinct  from  these 
methods  and  incorporates  light  transmission  though  an  entire  imaging  system  and  provides  calculations 
based  on  absorbance.  We  note  that  the  resolution  of  an  instrument  cannot  be  correctly  evaluated  unless 
the  pixel  size  is  appropriate  (smaller).  Even  if  two  objects  cannot  be  resolved  into  their  appropriate 
shapes  using  smaller  pixels,  the  detail  in  higher  pixel  density  images  is  higher.  Two  point  objects,  for 
example,  will  appear  as  dumbbells  or  ovals.  Hence,  even  for  systems  or  wavelengths  in  which  the  pixel 
size  is  smaller  than  required  for  correctly  sampling  the  spatial  frequencies  permitted  by  the  optics,  an  im¬ 
provement  in  image  quality  may  be  observed.  For  a  detailed  discussion,  we  refer  the  reader  to  appendix  A. 3. 

Data  from  the  two  optical  configurations  described  in  table  1  are  shown  in  Fig.  5.  Data  are  recorded  on 
prostate  tissue  that  is  prepared  using  methods  previously  reported23-25  and  images  are  obtained  using  the 
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Fig.  6.  Infrared  spectroscopic  imaging  data  from  breast  tissue  before  and  after  noise  re¬ 
duction.  (a)  An  absorbance  image  at  the  asymmetric  C-H  stretching  mode  before  noise 
reduction  and  the  corresponding  image  (b)  after  noise  reduction,  (c)  Spectra  corresponding 
to  recorded  data  and  data  after  noise  reduction  from  the  same  pixel. 


absorbance  of  the  asymmetric  C-H  stretching  vibrational  mode  (at  2962  cm-1).  The  top  row  demonstrates 
the  improvement  in  image  quality  attained  using  the  higher  NA  lens.  There  are  relatively  higher  spatial 
frequencies  present  in  the  data  acquired  with  the  higher  NA  lens,  signifying  that  there  is  information  loss 
in  increasing  the  pixel  size  from  1.115  fim  to  5.5  /im.  There  is  a  noticeable  set  of  high  spatial  frequencies 
present  in  the  data  acquired  at  the  higher  NA  at  the  shorter  wavelength.  Even  for  longer  wavelengths, 
the  information  content  difference  is  relatively  small,  but  not  zero.  To  capture  the  highest  image  quality 
across  the  spectrum,  hence,  the  pixel  size  should  be  calculated  based  on  the  highest  wavenumber  measured 
in  the  experiment.  This  statement  is  both  rigorously  and  intuitively  correct.  From  a  practical  perspective, 
however,  a  small  pixel  size  calculated  at  the  highest  wavenumber  also  reduces  the  throughput  significantly 
across  the  rest  of  the  spectrum.  To  address  this  trade-off  in  an  optimal  manner,  we  recommend  that  the 
highest  wavenumber  at  which  pixel  size  should  be  calculated  should  depend  on  the  typical  experiments  to  be 
performed  by  the  imaging  system.  This  small  adjustment  from  the  correct  sampling  at  the  highest  recorded 
wavenumber  to  that  at  the  highest  usable  wavenumber  for  image  generation  leads  to  the  concept  of  an 
optimal  pixel  size.  For  most  studies  in  the  mid-IR,  the  4000  —  400  cm-1  spectral  region  is  most  interesting. 
For  He-Ne  laser  reference-based  systems,  an  undersampling  of  the  reference  signal  by  a  factor  of  4  typically 
implies  that  the  allowable  free  scanning  spectral  range  is  3950  —  0  cm-1.  An  optimal  pixel  size  at  the  high 
end  of  this  range  is  0.974  fim.  In  most  experiments,  especially  with  biological  systems,  vibrational  modes 
at  this  high  limit  are  very  rarely  encountered.  A  more  practical  high- wave  number  region  is  3400  cm-1, 
which  is  in  the  vicinity  of  the  absorption  peak  of  O-H  and  N-H  stretching-associated  vibrational  modes. 
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Therefore,  we  calculate  the  optimal  pixel  size  at  this  wavenumber  and  the  predicted  pixel  size  is  1.13  tim 
on  a  side.  Since  we  sought  to  implement  the  concept  of  optimal  pixel  sampling  on  a  commercial  imaging 
system  without  extensive  hardware  modifications,  the  measured  pixel  size  of  1.115  fim  can  be  deemed 
acceptably  close.  Our  suggested  optimal  pixel  size  is  approximately  4-fold  larger  than  a  similar  setup  using 
the  synchrotron.  It  is  notable  that  intensity  considerations  are  secondary  for  a  synchrotron  source-based 
system  due  to  the  exceptional  flux  and  a  pixel  size  of  0.54  micrometer  on  a  side  was  used.  Relaxing  the 
very  strict  condition  with  a  more  practical  calculation  here,  we  maximize  the  spectral  quality  when  using 
a  globar  source  without  comprising  on  the  image  quality  in  any  appreciable  manner.  The  image  quality 
presented  here,  is  likely  of  the  highest  quality  that  will  be  observed  in  commonly  analyzed  biomedical, 
materials  or  forensic  samples  regardless  of  the  source. 

Though  image  quality  is  improved  with  a  smaller  pixel  size,  there  is  a  corresponding  decrease  in 
throughput  if  the  same  source  and  fore-optics  are  employed.  The  approximate  25-fold  reduction  in  pixel 
area  between  the  two  configurations  here  implies  that  acquisition  time  would  need  to  be  increased  625-fold 
(other  factors  being  constant)  if  the  data  quality  is  to  be  recovered  by  signal  averaging. 26,27  Some  of  the 
loss  is  mitigated  by  increased  throughput,  as  proposed  using  a  synchrotron, 13  or  by  increased  integration 
time  of  the  FPA,  as  in  the  experiments  conducted  here.  We  observed  a  6-8  fold  decrease  in  recorded  signal 
when  using  the  higher  NA  lens,  compared  to  the  lower  NA  system.  Hence,  the  need  to  signal  average  is 
not  especially  drastic.  Further,  there  are  other  methods  to  increase  the  SNR.  We  have  previously  proposed 
computational  noise  reduction28,29  to  obtain  a  significant  gain  in  SNR  without  the  corresponding  increase 
in  data  collection  time.  The  utility  of  this  idea  is  presented  in  Fig.  6.  A  significant  increase  in  spectral  SNR 
(Fig.  6(c))  without  observable  loss  in  image  quality  (Fig.  6(a)  and  (b))  can  be  observed.  Thus,  a  desktop 
high- definition  IR  imaging  system  is  indeed  feasible  and  can  provide  both  excellent  spectra  and  spatial 
image  quality. 

5  Conclusions 

A  complete  theoretical  understanding  of  the  image  formation  in  an  IR  microscope  is  provided  using  a  rigorous 
theoretical  model.  The  model  was  used  to  predict  the  optimal  pixel  size  at  the  sample  plane  that  would 
provide  the  highest  image  quality.  Simulations  demonstrated  that  the  effects  of  the  higher  NA  systems 
arose  from  an  increased  acceptance  of  angular  frequencies  and  resulted  in  higher  resolution  images  whereas 
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optimal  pixel  sampling  demonstrated  dramatic  improvements  in  image  detail  for  a  specific  NA.  The  results  of 
simulations  were  validated  using  measurements  on  two  different  configurations.  A  table-top,  high-definition 
FT-IR  spectroscopic  imaging  system  was  demonstrated  by  minimally  modifying  a  commercial  system.  The 
resulting  data  using  this  high-definition  system  can  be  of  high  spatial  and  spectral  quality  using  conventional 
signal  averaging  and/or  emerging  signal  processing  methods.  The  development  of  a  theoretical  understanding 
and  its  application  to  design  of  microscopes  and  acquisition  of  high-definition  data  should  spur  improved 
applications  in  many  fields  where  IR  imaging  is  applied. 
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Fig.  7.  A  Cassegrain  model.  A  point  source  placed  at  the  focus  results  in  a  plane  wave. 
The  outer  and  inner  numerical  aperture  limit  the  angles  that  enter  the  Cassegrain. 


A  Appendix 

A.l  Derivation  of  the  transfer  function  of  the  Cassegrain 

We  model  the  Cassegrain  as  a  focusing  system  with  two  concentric  apertures,  one  outer  aperture  with 
NAout  =  sin(#ont)  and  a  second  inner  circular  obstruction  with  NA^n  =  sin(^n)  (see  Fig.  7).  The  focal 
length,  which  is  the  distance  between  the  Cassegrain  and  the  focal  point,  is  assumed  to  be  known. 

Consider  an  object  placed  a  distance  d\  from  the  Cassegrain  as  shown  in  Fig.  8.  Let  the  field  in  a 
plane  perpendicular  to  the  principal  axis  at  the  object  be  \Uo).  The  field,  |J7i),  just  before  the  Cassegrain 
is  given  by 

\Ui)  =  Kdl  \U0) .  (21) 

The  field  immediately  after  the  Cassegrain,  | C/2),  can  be  expressed  in  terms  of  \Ui)  and  an  operator,  Gc, 
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\ 


Fig.  8.  An  object  placed  at  distance  d\  from  the  Cassegrain  produces  an  image  at  position 
^2  on  the  other  side. 


Fig.  9.  The  Cassegrain  is  flipped  and  the  object  distance  d\  is  larger  than  the  image  distance 
^2  on  the  other  side.  Equivalently,  a  larger  range  of  angles  are  accepted  by  the  Cassegrain 
from  the  object  side 
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corresponding  to  the  transmission  function  of  the  Cassegrain  as 


\U2) 

G  c 


=  Gc  \Ui) 


! 


d2r  |r)  Qc(r)exP 


.7T  V  9 
-i—rA 

Lf 


(r| 


(22) 

(23) 


where  Lf  is  the  focal  length  of  the  Cassegrain.  Note  that  the  quadratic  phase  function  in  Eq.  23  is  present 
because  we  have  used  a  thin  lens  like  approximation  for  the  focusing  system.  Qc{T)  =  Q^(r)  is  defined  in 
Eq.  24  and  R\  and  R2  are  the  inner  and  outer  aperture  radii. 


Qa(  r) 


{1  a  <  r  <  b 
0  else 


The  field  | C/3)  in  the  image  plane  is 


I  u3)  =  Kd2\  U2). 


(24) 


(25) 


The  term  fz( f)  in  ~Kd  can  be  approximated  as  fz( f)  «  v  ^1  —  ^2^.  After  substituting  Eq.  22  and  Eq.  21  in 
Eq.  25  and  simplifying  we  get 


I U3)  =  Kd2GcKdl\U0) 

J!  |m\  0^2/ 
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(26) 


This  is  the  relation  between  the  image  and  the  object  for  a  Cassegrain.  In  most  FT-IR  imaging  systems, 
we  can  make  a  few  approximations  and  simplify  Eq.  26.  The  object  size  (illuminated  area  of  the  sample) 


is  typically  much  smaller  than  the  size  of  the  Cassegrain.  This  means  that  exp 
size  is  also  much  smaller  than  the  Cassegrain  dimensions;  therefore,  exp 


//2 

1  2 d2  ' 


•  27 tv  ^/2 


1.  The  image 


1.  We  know  that30  for 


a  thin  lens  (or,  more  generally,  for  a  system  which  focuses  light  like  a  thin  lens),  the  position  of  the  image 


and  object  are  related  by  ^  .  Therefore  exp 


i“(A  +  5-r7)te2+!'2>]=1- 
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With  these  approximations,  Eq.  26  reduces  to 


| U3)  =  J d2rd2r'd2r"  |r)  Q^(r')  (r'\U0)exp 


27TZ/  . 

.2ttu  ,, 

i——r  •  r 
di 

exp 

1 - r  •  r 

d2 

(27) 


Defining  magnification  M  =  ^  and  evaluating  this  integral  after  projecting  on  the  Fourier  basis  yields 


ms)  =  QrI 


(28) 


It  may  be  noted  that  the  system  is  linear,  but  not  shift  invariant  (because  of  the  magnification  term).  With 
this  caveat,  we  can  think  of  Q ^  (— ^)  as  the  Cassegrain  transfer  function.  In  our  model,  the  exact  values 
of  focal  lengths  are  not  required  provided  that  the  system  is  in  focus.  Only  the  outer  and  inner  numerical 
apertures  of  the  Cassegrains  and  the  magnification  factors  are  required. 


For  a  Cassegrain  with  an  orientation  as  in  Fig.  8,  typically  the  image  is  close  to  the  focus.  Thus, 
we  can  make  approximations  that  d 2  ~  Lf  ( Lf  is  focal  length),  NAont  ~  and  NA^n  ~  Also  note 
that  Q  is  an  even  function.  This  gives 

(f\U3)  =  Q?NA°;‘(f)  (~Mf\Uo)  ■  (29) 

Note  that  propagation  in  free  space  is  a  spatial  frequency  bandpass  Qg(f)  and  spatial  frequencies  beyond  v 
do  not  reach  the  detector.  For  a  Cassegrain  with  an  orientation  as  in  Fig.  9,  the  object  is  close  to  the  focus 
and  so  we  can  make  the  approximations  that  di  «  L/,  NAont  «  ^  and  NA^n  «  This  gives 


<f|%>  =  QmA°T  (Mf)  (—Mi\Uo) 

(-f/M\u3)  =  e^r(f)<f  1^0). 


(30) 

(31) 
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A. 2  List  of  Symbols 


Symbol 

Description 

V 

Wavenumber 

NA 

Numerical  Aperture 

r  =  (x,y) 

Coordinate  space  variables 

f  =  (/*,/„) 

Transverse  spatial  frequency  variables 

fz 

Longitudinal  spatial  frequency  variable 

1  uZi) 

Field  in  a  plane  at  z  =  Zi 

1 S) 

Sample 

Qci 

Aperture  function  of  Cassegrain  Ci 

MCi 

Magnification  of  Cassegrain  Ci 

Kz 

Propagation  operator 

Idi,d2 

Operator  for  the  interferometer 

HCi 

Transfer  function  (operator)  for  a  Cassegrain  Ci  in  focus 

GCi 

Transmission  function  (operator)  of  Cassegrain  Ci 

Table  2.  A  list  of  symbols  used  along  with  their  description 


A. 3  Point  Spread  Functions  and  Resolution 

Point  spread  functions  for  the  two  configurations  in  table  1  are  presented  in  Fig.  10.  Simulation  data  for 
two  point  objects  separated  by  four  different  distances  is  presented  in  Fig.  11  and  Fig.  12.  Data  in  Fig.  11 
corresponds  to  the  high  NA  configuration  in  table  1  and  Fig.  12  corresponds  to  the  lower  NA  configuration. 
All  data  is  at  3950  cm-1.  Axes  in  all  the  images  are  scaled  to  the  effective  size  of  the  detector-intensity- 
image  at  the  sample.  From  these  images  it  is  evident  that  two  points  separated  by  1  /im  cannot  be  resolved 
using  either  NA c2out  =  0-65  or  NA c2out  =  0-5,  whereas  they  can  be  resolved  in  both  configurations  at  a 
separation  of  4  fim.  However,  NA c2out  =  0-65  can  resolve  two  points  separated  by  2.4  fim  (which  is  the 
resolution  according  to  Rayleigh  criterion),  while  NA c2out  =  0.5  cannot.  Points  at  a  separation  of  3  /im 
can  be  just  resolved  using  NA c2out  =0.5.  Data  presented  here  illustrates  the  best  image  quality  obtainable 
from  systems  with  the  optical  configurations  in  table  1  without  discretization  (i.e  using  an  analog  detector). 
Discretized  data  obtained  at  the  optimal  pixel  size  (as  calculated  in  section  2.2)  have  all  the  information 
needed  to  construct  images  of  this  quality. 
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Fig.  10.  Point  spread  functions.  Data  on  the  left  corresponds  to  a  configuration  with 
NA  c2ont  =  0  .65  and  data  on  the  right  corresponds  to  a  configuration  with  NA c2out  =  0.5. 
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Fig.  11.  Detector  intensity  at  3950  cm-1  for  two  point  objects  separated  by  distances 
indicated  on  top.  Corresponding  profile  plots  at  y  =  0  are  presented  in  the  bottom  row. 
Cassegrain  C2  has  NA c2out  =  0.65.  Axes  are  scaled  to  the  effective  image  size  at  the  sample. 
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Fig.  12.  Detector  intensity  at  3950  cm-1  for  two  point  objects  separated  by  distances 
indicated  on  top.  Corresponding  profile  plots  at  y  =  0  are  presented  in  the  bottom  row. 
Cassegrain  C2  has  NA c2out  =0.5.  Axes  are  scaled  to  the  effective  image  size  at  the  sample. 
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